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ABSTRACT 


This report introduces the reader to radiation hydrodyna mics 


(RH) and discusses its application to fireballs in the 


atmosphere. After formulating the ba sic equations of RH I 


special attention is given to the radiative transfer problem. 


Several methods for solving the equations of transfer are 


touched upon but special emphasis is placed on the two 


strea m method with a frequency averaging procedure I which 


is specifically designed for use with finite zone sizes. A 


version of the FIREBALL code which utilizes this approach is 


described. The physics of fireballs is illustrated with the 


exa mple of a one kiloton detonation at sea level density and 


without interference from the ground .. Some remarks are made 


on scaling procedures for extending the results to higher yields 


and altitudes. Estimates are made of the validity of the models. 
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FOREWORD 


"Thermal radiation ,i is electromagnetic radiation emitted by matter in a 
state of thermal excitation. The energy density of such radiation in an en- 
closure at constant temperature is given by the well known Planck formula. 


The importance of thermal radiation in physical problems increases as the 


temperature is raised. At moderate temperatures (say I thousands of degrees 


Kelvin) its role is primarily one of transmitting energy; whereas at high 


temperatures (say I millions of degrees Kelvin) the energy density of the radi- 
ation field itself becomes important a s well. If thermal radiation must be 
considered explicitly in a problem I the radiative properties of the matter 


must be known. In the simplest order of approximation I it can be assumed 
that the matter is in thermodyna mic equilibrium "locally" (a condition called 
local thermodynamic equilibrium I or LTE) I and all of the necessary radiative 
properties can be defined I at least in principle. Of course whenever thermal 


radiation must be considered J the medium which contains it inevitably ha s 
pressure apd density gradients and the treatment requires the use of hydro- 


dyna mics. Hydrodyna mics with explicit consideration of thermal radiation is 
called "radiation hydrodynamics" . 


In the past twenty years or so I many radiation hydrodynamic problems 


involving air have been studied. In this work a great deal of effort ha s gone 


into calculations of the equilibrium properties of air. Both thermodynamic 
and radiative properties have been calculated. It has been generally believed 
that the basic theory is well enough understood that such calculations yield 
valid results I and the limited experimental checks which are possible seem to 
support this hypothesis. The advantage of having sets of tables which are 
entirely calculated is evident: the calculated quantities are self-consistent 


on the basis of some set of assumptions I and they can later be improved if 
calculational techniques are improved J or if better assumptions can be made. 


iii 


The origin of this set of books was in the desire of a number of persons 


interested in the radiation hydrodynamics of air to have a good source of 


reliable information on basic air properties. A series of books dealing with 


both theoretical and practical aspects was envisaged. As the series materialized, 


it was thought appropriate to devote the first three volumes to the equilibrium 


properties of air. They are: 


The Equilibrium Thermodynamic Properties of Air, 
by F. R. Gilmore 


The Radiative Properties of Heated Air, 
by B. H. Armstrong and R. W. Nicholls 


Tables of Radiative Properties of Air, 
by Lockheed Staff 


The first volume contains a set of tables along with a detailed discussion of the 


basic models and techniques used for their computation. Because of the size of 


the related radiative tables and text, two volumes were considered neces sary . 


The first contains the text, and the second the tables. It is hoped that these 


volumes will be widely useful, .but because of the emphasis on very high tempera- 


tures it is clear that they will be most attractive to those concerned with nuclear 


wea pons phenomenology, reentry vehicles, etc. 


Our understanding of kinetic phenomena, long known to be important and at 


present in a state of rapid growth, is not as easy to assess as are equilibrium 


properties. Severe limitations had to be placed on choice of material. The 


fourth volume is devoted to general aspects of this topic. It is: 


Excitation and Non Equilibrium Phenomena in Air, 
by Landshoff, et al. 


It provides material on the more important processes involved in the excitation 


of air, criteria for the validity of LTE and special radiative effects. 
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A discussion of radiation hydrodynamics was felt to be necessary and 


the fifth volume which deals with this topic is: 


Radiation Hydrodynamics of High Temperature Air, 
by Landshoff, Hillendahl, et al. 


It reviews the basic theory of radiation hydrodynamics and discusses the 


application to fireballs in the atmosphere. 


The choice of material for these la st two volumes wa s made with an eye 


to the needs of the principal users of the other three volumes. 


Most of the work on which these volumes are ba sed wa s supported by the 


United States Government through various agencies of the Defense Department 


and the Atomic Energy Commission. The actual preparation of the volumes 


was largely supported by the Defense Atomic Support Agency. 


We are indebted to many authors and organizations for assistance and we 


gratefully, acknowledge their cooperation. Weare particularly grateful to the 


RAND Corporation for permission to use works of F. R. Gilmore and H. L. Brode 


and to the IBM Corporation for permisSion to use some of the work of 


B. H. Armstrong. Most of the other authors are employed by the Lockheed 


Missiles and Space Company, in some cases as consultants. 


Finally I we would like to acknowledge the key role of Dr. R. E. Meyerott 


of LMSC in all of this effort I from the initial conception to its realization. 


We are particularly grateful to him for his constant advice and encouragement. 


Criticism and constructive suggestions are invited from all readers of 


these books. We understand that much remains to be done in this field, and 


we hope that the efforts represented by this work will be a stimulus to its 


development. 
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Preface 


This volume reviews the ba sic theory of radiation hydrodyna mics and 


discusses the application to fireballs in the atmosphere. The first chapter 


starts with a formulation of the basic equations and goes on to discuss 


schemes for translating these impossibly difficult equations into manageable 


computing procedures. As a companion to this chapter we have added 


Appendix A with a version of Hillendahl' s FIREBl\LL code I which runs without 


inputs of a classified nature. 


Chapter 2 deals with the .physics of fireballs. The main discussion is 


devoted to the description of a one kiloton detona tion at sea level. That 


section has nearly all been written by H. L. Brode of the RAND Corporation 


but a few passages have been added by the editor. One of these deals with 


opaque precursors to shocks whose Significance to the thermal output was 


noted by Hillendahl since the original version was written. The section on 


other yield and altitudes wa s also written by the editor. 


The summary chapter examines the reliability of the results and how 


this is affected by approximations I incomplete ba sic information and other 


deficiencies in the present state of the art. 


I would like to thank Dr. H. L. Brode for his contribution and the 


RAND Corporation for permis sion to include his work in this volume. 


R. K. M. Landshoff 
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Chapter 1. 
RADIATION HYDRODYNAMICS 


1. 1 
Introduction 


A nuclear detonation deposits a large amount of heat energy in the 


air around it. The heating phase is of relatively short duration since the 


energy arrives in the form of X-rays which come either directly from the 


surface of the exploding bomb or from the shock heated air in the immediate 


vicinity of that surface. 


Following the X-ray deposition the air approaches local thermodynamic 


equilibrium (LTE). The method of calculating the subsequent explosion 


history which is discussed in this chapter ignores this period where the 


air relaxes to LTE. 
Before we proceed we take a short look at the validity 


of that assumption. 
* 
The kinetics of relaxation processes has been discussed in Chapter 6, (4) . 


The relaxation time depends on the ambient air density and on the final 


temperature as shown in Fig. 6.1 (4). 


For a detonation at sea level practically all the energy deposited 


by X-rays gets stuck in a relatively small volume and raises the temperature 


to very high values. Under these conditions relaxation times are very 


short. For a detonation at a high altitude a sizeable fraction of the X-ray 


energy is deposited at large distances and produces a lesser temperature 


rise because of the inverse square drop of the flux density. The lower 


air density and the lower temperature both contribute to increase the 


relaxation time. 


* DASA-1917-4, from now on referred to as (4). 
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As an example, let us consider a detonation with an X-ray yield of 


10 20 ergs radiating like a blackbody with a temperature of 10 7 oK occurring 
-3 
at an altitude somewhat below 50 km where the air density is 10 
times 


less than at sea level. A crude estimate, using the asymptotic theory of 


Section 4.4, (4.) shows that about 10% of the X-ray energy is deposited at 


distances more than about 80 m where it produces temperatures less than 


12,OOOoK. In Fig. 6.1 (4) one reads off that the relaxa.tion time at that 


temperature and a density 
--.Q... = 10-3 
is 10- 6 sec. Within that sphere 
Po 
it takes less time and on the outside more time to relax the air to its 


equilibrium temperature. Thus 10% of the energy relaxes at a relatively 


slow rate and the assumption that one can ignore the relaxation period 


is not entirely justified in that case. 


The assumption of LTE is essential to the classical formulation of 


hydrodynamics. It means that the temperature is a well defined property 


of the fluid and that pressure and internal energy are known functions of 


density and temperature. Without LTE it would be much more difficult to 


formulate the conservation theorems for momentum and energy. 
* 
In the theory of radiative transfer (Chapter 2, (2) , which together 


with hydrodynamics accounts for the expansion of fireballs, LTE also is 


an assumption of major importance. Without it a quantitative prediction 


of the interaction between matter and radiation would be a hopelessly 


complicated problem. 


Despite the very important role played by radiative transport 


radiation does not as a rule account for a significant fraction of the energy 


density and the pressure within a fireball. Even for blackbody radiation, 


* DASA 1917-2, from now on referred to as (2). 
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which is not present unless the gas is opaque, this contribution is small 


unles s the temperature exceeds values like 25 eV. Temperatures of that 


magnitude are only maintained during the very early stages of fireball 


histories. In this period the fireball cools down by radiative expansion 


and this goes so fast that there is essentially no hydrodynamical motion. 


In formulating the hydrodynamic equations one can therefore ignore the 


energy density and pressure of radiation because by the time they get 


into the act they are indeed negligible. 


During the early period of fireball expansion twhere radiative transfer 


of energy is important) the shape generally appears to be almost spherical, 


at least at low altitudes where the size is small compared to the scale 


height of the atmosphere. Asymmetries which are hidden by the opaque 


outer layers may possibly occur due to instabilities at the bomb air interface, 


but we shall ignore these. Not much is known about such phenomena in 


any case and adding the complication of asymmetry would compromise the 


already complicated problem of treating radiation flow. In line with the 


current state-of-the-art we shall therefore discuss only spherically symmetric~l 


problems. 


1.2 
Basic equations of radiation hydrodynamics 


The differential equations for calculating fireball histories are the con- 


servation relations of ordinary hydrodynamics but with a rather complicated heating 


term in the energy equation. They can be written in either Eulerian or Lagrangian 


form. The two forms are chara cterized by a different choice of independent space 


variables. In the Eulerian system these are the coordinates in real space and in the 
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Lagrangian one they are coordinates which are tied to the particles 


of the fluid. In the Lagrangian system the coordinates in real space 


which describe the position of a specified particle are used as dependent 


variables. In the Eulerian system this is manifestly impossible and the 


motion is described in terms of the fluid velocity. 


The other dependent variables which characterize the thermodynamic 


state of the fluid are the same in the two systems and can be chosen from 


a set which includes the density 
p 
or its reciprocal the specific volume 


V 
I the pres sure 
P 
I the temperature 
T 
I the internal energy 
E 
I etc. 


It may be convenient to keep several of these variables in the equations 


but one must keep in mind that they are interrelated and that the thermodynamic 


state is specified by any two of them. 


The Lagrangian method is especially useful in problems with a high 


degree of symmetry where one needs only one coordinate to specify the 


position. Having restricted ourselves to spherically symmetrical problems 


we shall therefore adopt the Lagrangian approach. 


We define the Lagrangian radius 
r 
of a given particle as its radius 


at time zero, i.e. before it has started to move. The actual radius of the 


particle at any time is denoted by the capital letter R. The hydrodynamic 


problem is to find 
R(r, t) • 


If 
stands for the initial density the specific volume at any 


instant is 


v = ~ (R)2 
Po r 
(conservation of mass) 
(1.2-l) 
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Introducing the velocity 


u = aft 
~t 


the other conservation equations are 


~u = 


~t 
P 
I 


O 
(;R)2 ::. ..... 
~ (conservation of momentum) 


. 
VQ 
(conservation of energy) 


. 
where the rate of heating per unit volume 
Q 
still needs to be worked out. 


As it stands the energy equation has a serious defect because it does 


not allow for the entropy raise produced by a shock. To get around this 


(1.2-2) 


(1. 2-3) 


(1. 2-4) 


we adopt the method of Von Neumann and Richtmyer (l950) and add a pseudo- 


viscous pressure 


{ (t Z ;'\I) 0 (a UJa rl 
if Q.!:!.< 0 
~r 


q = 
(1 .2-5) 
if 
~u > 0 
~r 


to the regular pressure in Eqs. (1. 2-3) and (1. 2-4). The constant t 
has 


the dimensions of a lengthi it will be further specified when we go to 


finite difference equations. 
. 
The radiative heating rate 
Q 
at some point is the difference between 


absorbed and emitted power per unit volume 


(1.2-6) 
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The emitted power presents no problem because the blackbody intensity 


B 
is a known function of temperature. There is no angular dependence 
\I 


and the integral over frequency can be expres sed in terms of the Planck 


mean, Eq. (2,4-15), (2). One obtains: 


Q• 
= 4\l-p (J T4 
em 
(1.2-7) 


where 
(J 
is the Stefan Boltzmann constant. 


The absorbed power is much more difficult to evaluate because the 


calculation of the intensity 
I 
is a major task. To carry this out one should 
\I 
in principle solve the equation of transfer (Eq. (1.3-1» along every ray 


pas sing through the pOint in question and for all values of the frequency. 


One of the major difficulties of such a program arises from the fact that 


the optical properties of air in a large part of the relevant temperature range 


result mainly from transitions between molecular levels. The spectrum 


associated with the major band systems consists of an enormous number of 


lines and the absorption coefficient fluctuates from large values at the 


line centers to small ones between the lines. Because of these "windows II 


the radiation at some point generally comes from pOints along the ray which 


are an appreciable distance further back. This distance varies just as 


strongly with frequency as 
U 
I 
. \I 
itself and it is therefore not proper to use 


local averages of 
I-l 
I 


\I 
in a frequency interval containing, say I a few lines. 


Instead, it is in principle necessary to integrate the transport equation at 


a large enough number of frequencies within everyone of these intervals. 


This correct approach clearly demands an impossible amount of computational 
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effort which has to be avoided. There are two limiting situations where 


this can be readily done. The one situation arises fer a transparent medium 


where the optical depth 
I...l IL 
v 
(L being the size of the radiating region) 


is uniformly small compared to unity. In that case 
Iv 
is very much smaller 


than 


have 


B 
and one can neglect the absorption altogether. In that case we 
v 
Q = - Q 
whichweknowfromEq. (1.2-7). 
em 


In the opposite extreme of an opaque region for which 
I...l IL» 
1 
v 


one can simplify Eq. (1.2-6) directly. The heating rate can in that case 


be expressed in the farm* 


Q = 
.... 
..... 
'I • F 


The flux vector is given by Eq. (2.5-8) I (2) which we rewrite in the form 


F = 
4 
..... 
4 
3" "R a 'I T 


(l. 2-8) 


(1. 2-9) 


Tables and graphs of the Rosseland mean free path "Ror the related opacity 


** 
can be found in (3) 
I p. 12 I pp. 446 to 449 and pp. 622 to 625. 


The above method of treating radiative transfer was originally 


developed by Eddington nearly half a century ago. In its application it 


was however I limited to astrophysical problems where it was not coupled 


to hydrodynamics. An early discussion of the use of this so-called diffusion 


approximation to radiation hydrodynamics has been given by Magee and 


Hirschfelder (1953). The first calculations carried out with this method 


to appear in the open literature were presented by Marshak (1958). 


* 


** 


The operator V is defined in Eulerian space. In plane or spherical 
geometry it is well known how to express it in Lagrangian form. 
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1.3 
Average absorption coefficients 


In the temperature range where molecular transitions occur and 


where optical depths are neither uniformly s mall nor uniformly large one 


has to resort to approximation schemes. It is clearly necessary to apply 


some kind of frequency averaging which will do a fair amount of violence 


to the "correct approach" of solving the transfer equation for a few million 


values of the frequency. The basic mathematical problem is that one wants 


to average the product 
~v I 
~ which enters in Eq. (1.2-6) as well as in 
* 
the transport equation 


dlv = 
ds 


~ 
I (B 
- I ) 
v 
v 
v 


by equating the average of the product and the product of the averages, 1. e. 


(1.3-1) 


one wants to replace 
~ 
I I 
v 
v 
by 
~ 
I I 
v 
v 
and that is of course not correct. 


The quality of this approximation depends on the amount of fluctuation among 


the values of 
~ 
I 
V 
and 
I v 
that are being averaged and in a line spectrum 


this fluctuation may be quite severe. A number of averaging schemes have 


been proposed and are used in various computing programs. 


One scheme divides the spectrum into groups (10 to 100) whose widths 


are chosen fairly narrow at the low energy end and wider as the energy goes 


up. Within each interval a Rosseland type average is obtained. Such group 


averages have been used, for example, in the SPUTTER program of AWFL as 


reported in RTD-TDR-63-3l28 Vol. II and in a code developed by J. Zinn 


of LASL. 


* 
Strictly speaking I the left hand side of this equation should contain the 


oI 


additional term 1 T 
but because of the large value of the light velocity 
this time depende~Je is usually left out. We note further that light rays are· 
straight lines in Eulerian space and in this section we temporarily abandon 
the Lagrangian system. 
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A second method of averaging uses the average transmission function 


(Eq. (2. 6-12b), (2) 


Tr {J.! '\Is) = _1_ J 
~\li 


~vi 


-~' s 
e 
v 
d\l 


and the slab absorption coefficient related to it by Eq. (2.6-19) I (2) 


~\I (s) = - 
1 In Tr (~ , s) 
s 
v 


These averages are defined for slabs of thickness 
s 
in which the 


temperature and density of the air are uniform. The intervals 
~\li 
are 


much narrower than the groups of the first mentioned method. The spacing 


between intervals is 10 to 20 times as large as the interval size. The 


calculated averages depend smoothly on the frequency so that it seems 


reasonable to interpolate. The slab average is made to order for use in 


finite difference equations where the fluid is divided into zones. It has 


been used in a number of LMSC codes which will be discus sed later in this 


chapter. 


A variation of the group average procedure consists of subdividing 


the frequencies within a group into subgroups which are ordered according 


to the magnitude of the absorption coefficient rather than by frequency. 


The spread between the absorption coefficients within each subgroup is 


( 1.3-2) 


( 1.3-3) 


obvious ly less than between those in the entire group and subgroup averages 
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will therefore be more meaningful. To use subgroup averages we must also 


introduce individual intensities for each subgroup. Even the use of only 


two subgroups would improve the accuracy considerably. A short-cu t for 


the calculation of two subgroup absorption coefficients consists of fitting 


the average transmission function in the form {Eq. (2.6-46). (2). 


Tables of 


- 
~ IS 


Tr (~ IS) = 1 
{e 
1 
+ 
v 
2 


~I 
and 
1 
u 2 are given by Churchill et a!. (1963). We 


don It know of any code which has utilized this type of average. 


1. 4 
Solution of the equation of transfer 


Having obtained an average absorption coefficient which permits us 


to replace the average product 
~v Iv 
by the product of averages 
~v Iv 
* 
the transfer equation becomes 


dIv = 
ds 
Uv (s - f ) 
v 
v 


The formal integration of this equation along a ray is straightforward 


and leads to 


* 


hV ' - T ) 


(s:) e 
dr v 


o 


The absorption coefficient is still meant to include the correction factor 
for induced emission (Eq. (2.2-11a), (2) and the prime is left out for 
convenience of writing only. 
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(l .3-4) 


(1.4-1) 


(1.4-2) 


S 


TV = T ).) = f ii)·I) d. I 
= T (s ') 
T\) 
\) 


So 


The difficulty arises because one has to determine the value of this 


integral for all the rays through a given point to evaluate the rate of 


absorption of radiative energy at that pOint. This is required for carrying 


out the angular integration in Eq. (1.2-6). 


There are basically two approaches to this problem. One is the brute 


force approach to follow this program directly and to evaluate 
~ (5) 
along 


(1.4-3) 


a large number of rays. This approach has been used in the SPUTTER program 


with one tangential ray through the center of each zone. Fig. 1-1 shows 


how these rays are combined to obtain the various values of 
I\) 
at the 


center of zone 4. Of the 7 rays which are drawn 3 are redundant because 


of symmetry and one obtains 3 different values of 


in and 1 grazing the zone. 


I \) 
going out, 3 gOing 


In the other approach one defines certain moments, i.e. angular 


integrals of I\) 
which now depend only on the radius and not on the 


direction. To solve for these moments one integrates a system of coupled 


linear differential equations which are only approximately correct but which 


give exactly the right answer when one considers the limit where the 


diffusion approximation applies. Such schemes h ave been used Widely 


in astrophysics and are discussed in great detail by Chandrasekhar (1960) 


and Mustel (1958). Some of the more sophisticated schemes use a large 


number of moments but quite good results can be obtained by restricting 


that number to two and using only the outgoing and ingolng flux which 
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are defined as the integrals 


cos 8 > 0 


(cos 8 < 0) 


where 
8 
is the angle between the ray and the radial direction. 


We consider firs t 
the case of plane geometry where the medium 


is stratified in plane parallel layers. This geometry has been studied 


extensively by astrophysicists and applied to the radiative equilibrium 


in the outer regions of stars where it is indeed unnecessary to worry about 


the curvature. 
By treating the radial coordinate 
R 
as if it were a 


cartesian coordinate the angle e of a ray remains constant along the 


ray path. The optical path length between two surfaces is therefore 


simply the optical path length along the normal divided by I cos 8/ • 


We shall expres s this in terms of the optical depth conventionally 


defined by astrophysicists as the optical path length measured radially 


(1. 4-4) 


inward from the surface of a star (or in our case a fireball) Ii. e. the integral 


R s 


T (R) 
\) 
=J 
0: (R') dR ' 
\) 
(1.4-5) 


R 


To evaluate 
as given by Eq. (1.4- 2) for an ou tgoing ray we place 


far enough inside that the factor 
e 
is essentially zero; for an 


ingoing ray we start at the surface where 
1\)(so) = 0 
• 
The first term of 
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Eq. (1.4- 2) can therefore be left out in both cases. For the exponent 


under the integral one can write 


, 
\I 
\I 
I 


T' ' - T' 


T \I 
- 
T \I = - 
. cos 8 


and for ou tgoing and ingoing rays one obtains 


I = 
1 
\I 
1 cos 81 


r 
f 
B (R') ~ (R') 
\I 
\I 
I 
T" 
-T' I 
~ 
\lcos e\l dR' ; cos 8 > 0 


o 


r s f B)R'l ii)R'l 


r 


cos 8 < 0 


T' '- T I 
\I 
\I 


e 
cos e 
dR' 


Entering these expressions into Eq. (1.4-4) one obtains the outgoing and 


the ingoing flux 


where 


r 


F\I+ = 2TT f B)R') I1)R') E2 ~T\I' - 'T)) dR' 


o 


r s 
Fv_ = 2IT f Bv(R'l il)R'l E2 ~1'v' - 1'vI) dR' 


r 


CD 


E2(T) = f e-UT u- 2 du 


1 
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(1.4-6) 


(1.4-7a) 


(l.4-7b) 


(1. 4-8a) 


(1.4-8b) 


(1.4-9) 


A useful approximation is obtained if one replaces 
/ cos 8/ 
in 


Eq. (1.4-7) by an average 
c = cos e. Subs ti tu ting the approximate 


form of 
Iv±' 
into Eq. (1. 4-4) gives expres sions similar to those in 


Eqs. (1.4-8) but the exponential integral is now replaced by a simple 


exponential function, i.e. we have the approximation 


The average intensities calculated from the approximate fluxes 


satisfy the differential equations 


dIV± 
-- 
c 
dR 
+ II (B - I ) 
- 
i'-"v 
v 
vi. 


These average intensities are therefore identical with the intensities in 


the directions for which 
/ cos e / = c. The idea of the two stream model 


with intensities in a characteristic direction goes back to Schwarzschild 


and Schuster who suggested to use 
1 
c = -2 
is 
c = t which gives the correct net flux 


in the high opacity limit. 
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• A much better choice 


dB 
~ 
dR 


(1.4-10) 


(1.4-11) 


(1. 4-12) 


(1.4-13) 


In spherical geometry one has no s1mple r1gorous expressions for 


I 
like those given in Eqs. (1.4-7) from which to derive two stream 
\I 
equations. It seems I nevertheless I reasonable that one should be able 


to use equations which are essentially of the same character but with 


minor modifications to maintain conservation of energy. It is easy to 


see that this is achieved by the pair of equations 


With the definition given in Eq.( 1.4-11) one obtains the outgoing and 


incoming flux simply by multiplying the corresponding intensities 
1\1+ 


by a factor 
TT 
• 
From the total integrated net flux 


:1 = f 
(F + - F 
) d\l 
\I 
\1- 


one can finally obtain the heating rate 


Q= 


for use in the energy equation. 


1. S 
Finite difference equations 


It is obviously impossible to find exact analytical solutions to the 


equations of RH and one must be satisfied with approximate numerical 


(1.4-14) 


(l.4-1S) 


(1.4-16) 


solutions. To obtain these one replaces infinitesimal increments of dependent 


as well as independent variables by finite differences. Mathematically RH 


IS 


can be characterized as an initial value problem and the methods and problems 


arising in treating this by means of finite difference equations have been 


. thoroughly discussed by Richtmyer (1957). We shall review some general 


considerations and then turn to ques tions which are specifically relevant 


to our problem. 


In a finite difference scheme continuous variables are replaced by 


discrete ones but there are numerous possibilities for doing this. Thus 


one can regard the discrete values of a variable as representing either 


the values of the corresponding continuous variable at a set of discrete 


meshpoint or the average values between meshpoints. There are other 


variations but they are not needed in the following discussion. One can 


treat some variables in the first and others in the second manner. To 


indicate the actual choice one can use integral subscripts for variables 


defined at the meshpoints and half integral ones for those defined in the 


intervals. It is convenient and natural to let 
Ri 
and 
Ui 
represent 


the radius and the velocity of the particle at the meshpoint 
1 
and this 


leads almos t automatically to defining 
V i+ 1/2 
I Pi+ 1/2 
I £i+ 1/2 
I 


and 
T1+ 1/2 
as the averages of specific volume I pres sure I internal 


energy densi ty and temperature in the interval between the meshpoints 


i 
and 
i+ 1 • 


Another element of choice enters in the methods used for advancing 


variables in time. Almost all variables are defined at meshpoints in time 


which are indicated by integral superscripts. It may I however I be useful 


to define the velocity between meshpoints which can be indicated by half 


integral superscripts. With this definition and abbreviating the right hand 


side of Eq. (1.2-3) by 
a 
(for acceleration) that equation and Eq. (1.2-2) 
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lead to the integration procedure 


Having obtained 
Rn+ lone can then obtain 
Vn+ 1 
by differencing 


which follows from Eq. (1. 2-1). So far we have not bothered to look at 


alternate schemes because the procedures outlined above are very straight- 


forward and there seems to be no good reason for doing anything more 


elaborate. In the purely hydrodynamic case I i. e. if Q = a the energy 


equation (1.2-4) can also be integrated Vf~ry simply. Centering the 


difference equation at (n+ 1/2) leads to 


If 
E 
is expressed as a function of 
V 
and 
p 
this equation can be 


n+l 
solved for 
p 
I the one variable which is still unknown. Anticipating 


the problems which arise when one has radiative heating it is really more 


useful to expres s both E 
and 
p 
as functions of 
V 
and 
T 
and to 


solve Eq. (1.5-3) for 
Tn+1. Either way one has to solve for only one 


unknown at a time which causes no real difficulty even though it may have 


to be done by iteration. 


This situation changes drastically when the radiative heating rate 
Q 


becomes important. The heating term to be added on the right hand side of 


Eq. (1.5-3) should be centered at the level 
tn+ 1/2 
like the remaining 
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{ 1. 5-1} 


{l. 5- 2} 


(1.5-3) 


part of the equation but the temperature distribution from which it must 


be calculated is only known up to the time 
n 
t 
• There are two major 
. 
avenues of attack. One is to forget about centering 
0 
and use its 


value as calculated at tn. If this is done one can still solve explicitly 


for 
TnT 1 
and this is as the explicit method of integration. 


In the other attack one uses the properly centered heating rate 


1 (. n+ 1 
• n) 
"2 \0 
+ 0 
. This means that the equation which describes the heating 


in anyone zone depends on the values of 
Tn+ 1 
in all zones so that one 


has to solve a large number of equations (one per zone) simultaneous ly. 


This implicit method involves a considerable amount of algebraic labor. 


If centering was only required for accuracy it would not be worthwhile to 


go to all this trouble because one could increase the accuracy more easily 


by reducing at. What is really involved is the question of mathematical 


stability which we shall briefly discuss. 


It is physically clear that a fluid responds to any pressure or temperature 


disturbance by a motion or heat flow which counteracts the disturbance. In 


an integration by means of difference equations which uses too large time 


intervals it may happen that the disturbance is overcompensated so that an 


excess tur ns in one step into a deficit, in the next step again into an excess 


etc. If the magnitude of this alternating disturbance increases each time any 


small disturbance will eventually cause the solution to blow up. In principle 


one can cure such an instability by taking 
at 
small enough but this could 


seriously increase the running time of a problem. 


There are two cases where the s tabil1ty condition has been obtained 


analytically. The firs t arises when the dominant mode of energy transfer 


is of a hydrodynamic nature. The maximum 
a t 
in this case is found as 
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follows. One calculates for each zone the traversal time 
t:. tiT 1/2 = 


(Ri+ l-RrV s (i+ 1/2) 
of a signal traveling with the local sound speed. Going 


through all intervals one then finds the smallest, say 
t:.tmin • The time 


increment is then limited by the so-called Courant-condition 


where 
k 
is a numerical factor near unity which depends on the integration 


scheme. In the scheme where one uses the three equations at tre beginning 


of this section one has 
k = 1. 


When radiative heating dominates, the stability analysis has been 


carried out for the case where one can use the diffusion approximation. In 


the explicit scheme the limit for 
a t 
is proportional to 
~Ra R2 
which 


decreases together with the Rosseland mean absorption coefficient of the 


air. If the air is fairly transparent 
a t 
is limited to very small values 


and this makes an explicit calculation very costly in computer time. The 


implicit method does not have this trouble and is in fact unconditionally 


stable. On the other hand it is of course also time consuming to solve a 


large number of coupled equations simultaneously. One can attempt to 


approach the implicit solution by iteration. On the first go-around one 


can advance 
T 
by the explicit method. With the advanced temperature 


Oon+l 
Oon+l/2 __ 
distribution one can then work out 
, form the average 


1 I· n 
• n+ 1) 
n+ 1 
2" \0 + 0 
and reevaluate 
T 
• This procedure can be repeated 


several times and if it converges it will lead to a stable solution. The time 


step 
a t 
is now limited by the condition that the solution should converge. 


In contrast to the stability condition of the explicit method this limit of 
at 
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(1.5-4) 


is inversely proportional to - 


~R and independent of the zone size. The 


actual convergence criterion is almost equivalent with imposing a limit on 


the fractional energy change per time step within every zone. That form 


of the condition is easy to use and experience has shown that a fraction 


like one percent ensures the convergence. In a purely implicit procedure 


there is no such limitation on the magnitude of the time step. For the 


sake of accuracy one should also impose a limit on the fractional energy 


change per time step but it does not need to be as small. This limit can 


be allowed to vary from zone to zone to require greater accuracy in those 


zones where the changes make a significant contribution to the overall picture. 


It is obvious that the allowed time interval changes throughout the 


calculation. For reasons of economy one should always run fairly close to 


the maximum without, however increasing 
Ii t 
too abruptly. To change 
<5 t 


generally requires some interpolation (or extrapolation) and all programs 


nowadays have provisions for carrying the necessary changes out automatically. 


Although the preceding arguments were based on the diffusion approximation 


they apply equally in the more general case. It is true that one can not 


readily obtain analytic stability or convergence criteria but experience 


with numerical calculations indicates the same pattern. 


In addition to the various decisions described above one also has to make 


a choice on zone sizes. There are two parts to this decision relating to the 


total number of zones and to their relative sizes at different radii. Part one 


involves a compromise between conflicting requirements for accuracy and economy 


because it takes a large amount of computer time to use very many zones. 


This is amplified if the choice of 
15R 
also limits the time step as in 


hydrodynamic calculations where 
Ii t '""" 6 R 
and even more in the explicit 
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calculation of radiative transfer where 
2 
ot.-voR. 


Part two involves a judg ment as to where the really significant 


changes are taking place and it is of course at those regions where one 


should use the finest zoning. In the course of a calculation the location 


of significant changes moves so that one has to make provisions in the program 


to detect this and to react to it by rezoning. Furthermore the overall radius 


of the fireball changes during an average calculation by as much as 3 orders 


of magnitude so that rezoning is also necessary to keep the number of zones 


at a more or less constant level. 


The pseudo viscous pressure 
q 
introduced in Eq. (1.2-5) is a device 


for calculating the entropy rise behind the shock. Without the damping 


mechanism provided by 
q 
the changes induced behind the shock overshoot 


and produce lasting oscillations which physically do not belong there. A 


large value' of 
t 
will kill these improper oscillations most effectively 


but at the cost of making the transition region very wide which is also 


incorrect. Experience has shown that 
.(, = 26R 
will stop the fake 


oscillations reasonably fast without spreading the shock transition over 


more than about 4 zones. 


One can also express 
q 
in terms of 


Richtmyer suggests to use the formula 


with 


(p .(,)2 


q _ ---=o~_ 


- 
V 
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o v rather than 
ot 


oV < 
0 
at 


~ 
or 
and 


(1.5-5) 


(1.5-6) 


SO that the transition region covers the same number of zones near the center 


and further away from it. The numerical factor 
a 
should again be 


approximately 2. 


In differencing either Eq. (1.2-5) or (1.5-5) one is led to expressions 


at half integral times. To obtain the acceleration in Eq. (1. 5-1) it should 


be known at 
tn 
but to achieve that, one would have to use an implicit 


rou tine. In this case that is not worthwhile since the use of 
q 
is an 


artifice anyway and it is customary to have 
q 
lag half a time step behind. 


In the energy equation (1. 5- 3) 
q 
is automatically in step. 


The total energy which is obtaired by summing the kinetic and internal 


energy within the fireball and the energy carried away by radiation should 


always stay at a constant level. The internal energy should in principle 


contain a part due to radiation but as mentioned in section 1.1 this does 


not amount to much. A trivial point, but one which must nevertheless be 


kept in mind is, that one should only count the excess over the energy in 


the ambient unheated air; otherwise the nominal energy would grow with 


the volume of the fireball. 


It is important to keep track of any violations of energy conservation 


which may creep in through the use of finite difference schemes. Any 


program should therefore contains a routine for checking energy conservation. 


The pOint at which R. H. goes beyond standard methods comes with 


the calculation of radiative transfer. The various methods require the 


evaluation of certain space integrals before one can calculate the energy 


deposition in a specified zone. Because of the very strong temperature 


dependence of the integrands these integrals depend critically on the 


radial dependence of the temperature. The common method of approximating 
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this dependence by assuming constant values of the temperature within the 


zones may lead to serious errors. An attempt to correct these has been 


made by Hlllendahl (1964). 


We shall present the analysis for the plane case which is formally 


easier. The transition to spherical geometry can be made later and 


requires only minor changes which are rather obvious. The starting pOint 


is Eq. (1.4-12) but before one has carried out any frequency averaging. 


Thus the line character is still preserved and 
Ilv 
and 
Iv+ 
are 


rapidly changing functions of frequency. Integrating Eq. (1.4-12) 


across the zones which are separated by the interface at 
Ri 
one finds 


for the outgOing and ingoing stream (represented by the upper and lower 


sign 


-].6 
-1 
2 \i i+- 
, 
2 
3 
I 
= I + i+-leT 
- 
v+,i 
v_, 
- 
2 


r v,i 


1 


I 
B v 


3 
1 
-- t. 
2 v,i 
..... 1 
e 
dT v 


I 


where 
B \i 
and 
are taken at the same pOint 
R 
, and 


1 
6 v,i 
t.v ,i + 1/2 = 


To carry out the integral we use the first two terms of the power expansion 
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(1.5-7) 


(l.5-8) 


(1.5-9) 


and obtain 


I v+ ,i 


-{- 1 (~) f -(1 + d. 6 
- 
) 
- ~ 6v , i -+ 1/2] 
- 
3 
c1T 
2 
v Ii -{- 1/2 e 
v 


(1.5-10) 


Integrating over frequency we obtain formally 


Z-{- 
i 
+ (Bi 1 - A+ 1) -{- i(d~) W+ i 
-' 
-' 
dT-' 
i 


(1.5-11) 


where 


(1.5-12) 


Z-{- 
i = J 
_I 


I v + , i 
1+ i 
_I 


3 
- 
/ 
-2'6,i+12 


e 
dv 
(1.5-13) 


(1.5-14) 


(1.5-15) 


and where 
Bi 
is the integrated intensity 


B = f 
B dv = Q.. T4 
v 
IT 
( 1.5-16) 
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at the pOint 
Ri • No limits of integration have so far been specified 


and one is free to divide the spectrum into any set of frequency intervals. 


As a first approximation H1l1endahl used the entire spectrum without 


subdividing it. 


Before Eq. (1.5-11) can become operational one has to define 


the average optical depth T which enters in the derivative 
and 


one has to face the difficulty of an unknown ratio 
I /1 
entering into 
v 


the definition of Z. 


The procedure devised by Hillendahl for obtaining an average for 


T is specifically intended for use with finite zone sizes. The 


prescription is designed to keep the emissivity of a zone of constant 


density and temperature unchanged if one replaces the frequency 


dependent optical depth 
t:. Tv = I-lvt:.R 
by its average 
t:. T • Thus, 


i.e., by making this substitution in the exponent of Eq. (1.5-14) one 


is led to 


A= J ~ e 


3 
-- ~ t:.R 
2 v 
dv = e 


3 
,...., 
-- t:. 'T 
2 


B 
We note that the factor ....Y..J... 
in Eq. (1.5-14) is taken at the edge of 
Bi 
the zone. Since this ratio varies only very s lowly with 
T 
we will 


take it at the center of the zone instead, so that there is only one 


emissivity 
A per zone and not different ones for the ingoing and out- 


going ray. In the above integral for 
A one can clearly replace the 


rapidly varying exponential by a smooth one in which one uses the slab 


absorption coefficient 
~v as defined in Eqs. (1. 3-2) and (1. 3-3). In 
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(1.5-17) 


the new expres sion 


3 - 
R 
-"2 \1 v !::,. 


e 
dv 


the zone with 
!::,.R 
enters not only as a factor in the exponent but also 


as one of the variables in 
0. = ~ (p, T, !::,.R) • Various methods for 
v 
v 


calculating 
~ ,which apply when the dominant absorption is molecular I 
v 


atomic or free-free, have been described in (2). 
The results are tabu- 


lated in and have been used in the above integral to obtain 
A{p IT ,!::,.R) • 


It is convenient to express this in terms of a mean absorption coefficient 


0H (p, T ,LlR) = - ~ tn A/!::,.R 
and to write: 


A= e 


_1 f,~ 
LlR) 
2 \f-LH 


Depending on 
!::"R 
as well as on 
p 
and 
T I 
~H differs from the 


Rosseland mean (il R) 
and the Planck mean C!1p) 
which depend only on 


p 
and T. 


For the function W, which should in principle be calculated from 


Eq. (1. 5-15) we use an approximation and set 


( 
3 
) 
- ~ ~H LlR 


W = 1 - 
1 + "2 ~H !J. R e 


which looks reasonable and leads to the correct energy deposition when 


~H!::"R is large enough that one can use the diffusion approximation. 


As in the case of 
A 
we are using only one 
W 
per zone. 
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(1.5-1S) 


(1.5-19) 


(1,5-20) 


The coefficients 
Z... i 
as defined by Eq. (1.5-13) depend on 
-' 
the unknown spectral distribution 
1\)+ ,i/1+,i at the pOint 
i • 


Since the difference equations (1.5-11) apply only to the integrated 


intensities 
1-+ i 
their solution does not give us any direct information 
-' 
about the spectral distribution and we must rely on educated guesses 


for the latter. The basic clue which we follow is that the radiation 


at some point 
i 
comes by and large from a zone (the radiating zone) 


which lies an optical depth unity behind that pOint in the direction 


where the stream comes from. The distribution has therefore the distri- 


bution of a blackbody source at the temperature 
TR 
of the radiating 


zone but modified by selective absorption in the intermediate zones. 


When we apply this model we dis tinguish 3 typical situations 


for the ingoing and 3 for the outgoing stream. This comes from the 


peculiar te'mperature dependence of mean absorption coefficients. For 


all these means I whether we talk of - 
- 


~R ' ~p or 
~H ,one can 


distinguish a central temperature range where - 
~ is large and the 


low and high 
T 
ranges where it drops to very low values. 


The temperature profile of a fireball is typically a monotonically 


decreasing curve. While the central temperature 1s still large this 


profile looks like the sketch in Fig. 1-2 with a central section where 


~ is small so that radiative transfer maintains a nearly constant 


temperature. Beyond t hat plateau comes a more opaque region with a 


relatively large temperature gradient and at the pOint where the 


temperature has dropped to where the air is again transparent the profile 


becomes again more level. Superimposed on this one finds usually 
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some structure due to shocks or other disturbances but this does not 


alter the main conclusion that there are 3 distinct regions. 


In the opaque intermediate region the spectral distribu tion 
I /1 
\J 
can be identified with that of a blackbody at the local temperature so 


that one is led to 
Z+ = Z_ = A • In the interior region the radiation 


comes mainly from its boundary where it is in contact with the opaque 


region. Since the temperature profile is quite level one does not commit 


a significant error by identifying 
I /1 
again with the local 
\J 
B /B 
\J 
which varies much less with 
T 
than either 
B 
or 
B itself. As 
\J 
in the previous case we therefore use the approximation 
Z = Z = A • 
+ 
- 


Only in the outer section do we have to make a more careful choice 


of 
I /1 
and only for the outgoing stream. The ingoing stream carries 
\J 
essentially no energy and it doesn't matter much what one does. The 


simplest choice is again to set 
Z_ = A • 


The pOint where it really counts that our model should adequately 


represent the true physical nature of radiative transport comes when we 


consider the outgoing stream as it emerges from the opaque region. 


The absorption to which this stream is subjected is largely due to 


molecules.. In calculating which parts of the spectrum are and which 


are not transmi tted one is greatly helped by the character of the energy 


dependence of 
~ • Fig. 1-3, which is a typical example taken from 


SACHA type calculations shows that 
~ is a very rapidly rising 


function of frequency.. From ,this graph we find by inspection that a 


zone of abou t 10 m thickness would transmit practically no photons 


above 5 eV and practically all photons below 3.8 eV.. Approximately 


one can assume that there is a critical photon energy 


28 


h\J 


C 
in the 


vicinity of 4.4 eV at which the transmitted flux is sharply cut off. For 
B 
a stream which starts out as a blackbody spectrum 
~ = bv (TR) 


one finds that the transmitted fraction of the energy is a known integral 


Vc 


L (TR, hv cJ = J 
bv (TRJ dv 


o 


The procedure for determining 
Z+ 
for any zone 
i - 1/2 
outside 


the opaque region starts out with finding the radiating zone which belongs 


to it and whose temperature has been designated as 
TR • We assume 


that the model spectrum 
I 
starts out there as a blackbody spectrum 
v 


Bv (TR) • Any of the zones through which it passes will not transmit any 


radia tion above its cu t-off energy 
hv c 
and the model spectrum which 


finally enters into zone 
i - 1/2 
remains 
I = B (TR) 
up to the 
v 
v 


lowest cut-off energy 
h\)min 
encountered by the stream but is reduced 


to 
I\) = 0 
above 
h\)min. If zone 
i - 1/2 
has a lower cut-off 


h\) c, i- 1/2 
we set therefore 


= L(TR h\) c, i- 1/2) 
L(TR h\)min) 


If the cut-off energy is equal to or larger than 
hV min 
the model 


spectrum would lead to 
Z+,i = 1 
which can It be right and indicates 


that the sharp cut-off approximation is too crude to fit this case. An 


upper limit for 
Z can be obtained by setting it equal to the local 
A 


since the bulk of the true spectral distribution lies at somewhat higher 


energies than the local blackbody spectrum. Thus the true 
I 
will 
v 


suffer somewhat more absorption than 
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B \) 
which leads to the inequality 


(1.5- 21) 


(1.5-22) 


Z < A • 


The modifications which are necessitated by going to spherical 


symmetry can be written down without difficulty. One only has to note 


that Eq. (1.4-14) differs from Eq. (1.4-12) by the factor 
R2 
which 


multiplies both 
I 
and 
B. Clearly this factor must enter when one 


modifies the corresponding set of Eqs. (1. 5-11). As we write down the 


modified set we incorporate the result that the coefficients 
A 
and 
W 


depend only on the zone and not on the direction of the stream, and obtain: 


The 
A 
and Ware as before given by Eqs. (1.5-19) and (1.5-20) and 


Z+ i 
is nearly always equal to 
Ai + 1/2 
except in a few zones just 
-' 
outside the opaque region where one should use Eq. (1. 5-22) •. 


In the two stream model the fluxes differ from the intensities by a 


factor 
TT 
as shown in Eq. (1.4-11). The relation carries over when 


one performs the frequency radiation SO that 


= 


Having determined the intensities by solving the set of Eqs. (1.5-23) one 


is therefore ready to evaluate 
.... 


- 
'V • F 
which according to Eq. (1.2-8) 


gives us the radiative heating rate. Thus one obtains: 
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(1.5-23) 


(1.5- 24) 


(1.5-25) 


which is the final equation of the two stream method. Before leaving it 


let us take a somewhat closer look how this equation handles a zone in 


the opaque region. In any region where one replaces 
Z+ 
and 
Z_ 
by 


A 
one finds that Eq. (1.5- 23) reduces to: 


from which one obtains in turn 


(1.5-26) 


(1. 5-27) 


In an opaque region this simplifies still further because the factor 


-Iii 6R 
2 t-"H 
e 


in Eq. (1.5-20) becomes negligible compared to unity. One can therefore 


set W = 1 
and obtain: 


and since 
B = (J T4 
this is clearly equivalent to Eq. (1.2-9) i. e., to 
n 


the basic equation of the diffusion approximation. This is of course no 


2 
surprise because we picked 
cos e = '3 
precisely in order to achieve this 


equivalence. 


In an opaque region, the diffusion approximations contains all 


the physics needed for the calculation of radiative energy transfer and 


it 1s superior to other methods as far as speed and possibly accuracy of 


calculation are concerned. Hlllendahl ' s formulation of the two stream 
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(1.5-28) 


method au tomatically leads to this procedure'. More elaborate methods 


such as the multiple ray technique do not, but it is of course possible 


to switch to a diffusion theory calculation when one considers an opaque 


region. Thi s is indeed done in the SPUTTER program. 


In the form outlined in this section the two stream method is not 


applicable at high altitudes where the ambient air has a density less 


/ 
-4 
than about 
p Po = 10 
• At such densities the air becomes transparent 


in the spectral region where 
B 
has its maximum and the fireball has 
\J 
no opaque region. There is still a significant amount of radiation at 


frequencies above and below this window but one has to devise new 


methods for dealing with this problem. At these altitudes the mathematical 


difficulties are further aggravated by the non-spherical energy deposi tion 


which takes place when the mean free path of x-rays get large compared to 


the atmospheric scale height. Eventually, say at about 
-6 
p/po~ 3 x 10 


the air becomes transparent in all parts of the spectrum and thermal 


radiation is no longer a significant factor. 
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FIG. 1-1 
RAYS CONVERGING UPON ZONE 4 WHICH ARE USED TO 
COMPUTE 
I 
AS FUNCTION OF ANGLE 
\J 
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Flli. 1-2 
SKETCH OF TEMPERATURE PROFILE IN A FIREBALL 
INDICATING OPTICAL PROPERTY OF THE AIR 
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FIG. 1-3 
TYPICAL ENERGY DEPENDENCE OF ABSORPTION COEFFICIENT 
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Chapter 2. 
THE PHYSICS OF FIREBALLS 


2 • 1 
In trodu ction 


A nuclear explosion in the atmosphere creates a fireball whose 


development depends in large measure on the physics of hot air. All 


of the previously discussed properties of hot air and all of the mechanisms 


for energy transport developed in previous chapters are a part of nuclear 


fireball physics. However I these energy transforming and transporting 


relations and the detailed knowledge of the properties of air find 


considerably wider application. They have or can contribute to the study 


of stellar dynamics I the nature of stellar atmospheres I the radiation 


from various astrophysical sources, and they can aid in the study of 


hypervelocity flight, upper atmosphere physics, aurora, and other atomic 


and molecular physics problems which involve high temperatures. 


It is certainly the case that the information presented in these 


previous chapters makes the conditions created in a nuclear explosion 


more understandable. Some knowledge of air heating mechanisms, of 


air excitation, of radiation transport, and of hydrodynamics, of absorption 


properties, and of the thermodynamics of air is necessary before a full 


description of a nuclear explosion can become more than heuristic. 


Much of the present knowledge about fireballs has been gleaned 


from test observations, but by far the greatest detail has come from 


numerical computer calculations, as have the quantitative estimates of 


fireball interior dynamics which appear in this chapter. Calculations of 


widely varying detail and sophis tication now abound, and it is not the 


intention in this chapter to review such results or analyze computing 


methods. Most current calculations rely for their measure of·success on 


the extent to which the physical concepts and properties covered in the 
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preceding chapters have been taken into account in a mathematical model. 


The principal objective of this chapter is to outline the physical features 


of nuclear fireballs and their thermal radiations, stressing where possible 


those factors which are most general and which provide the best under- 


standing on which to base predictions and extrapolations. The approach 


adopted is to begin by considering a small yield explosion (1 kiloton) 


at sea level and to describe the sequence of events which occur un- 


encumbered with interactions from the earth I s surface or inhomogeneous 


environments. This development will then be extended to higher yields 


and altitudes. There will be no attempt at completeness and no great 


concern for quantitative rigor I but it is intended to display as much as 


pos sible the current understanding of the physics of nuclear fireballs. 


2.2 
One kiloton at sea level 


A one kiloton explosion in a sea level atmosphere provides an 


appropriate example for an initial examination of the sequence of events 


that constitute a fireball history. The now familiar usage of kilotonage 


and megatonage refers to the total energy release in a nuclear explosion 


with the usual metric prefixes for a thousand or a million and with the 


understanding that a ton of high explosive - TNT - releases 10 9 calories 


of effective energy I i.e., one gram of TNT is taken as equivalent to 


one kilocalorie or 4.185 x 10 10 ergs. 


For any nuclear explosion the sequence of events is remarkably 


complex. In following its development for this one kiloton sea level 


explosion, the reader may appreciate that the present understanding, 
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although not complete, has become quite detailed and much of it has 


grown directly from the material reviewed in this series of volumes. 


The nuclear energy is released in an extremely short time - a 


small fraction of a microsecond - and always in a small mass and 


volume. It is the properties of this small mass, constituting the 


weapon itself and its carrier that determine the early source of energy 


for the fireball, and some of these properties may influence the 


character of the later thermal radiation. Everything starts in this 


nuclear source and all of the initial radiations - gamma rays, neutrons, 


and x-rays - 
are generated by it. However, the air or other immediately 


surrounding material absorbs almost everything emitted within a few 


hundred meters and the nature of the observable fireball is largely 


determined by the properties ot'this surrounding air. For our exa mple 


of one kiloton in a sea level atmosphere, the air within a few feet of 


the weapon stops nearly all of the x-rays, and the prompt gamma rays and 


the neutrons have removal mean free paths of about 400 and 240 meters, 


respectively. These rapid absorptions make knowledge of specific 


details of the nuclear device largely unneces sary in describing the 


fireball phenomena. Consequently, we shall be able to proceed without 


reference to classified aspects of nuclear weapons and yet without 


significantly truncating our description of the fireball and its thermal 


radiations. 


The fraction of the energy which may be radiated out of the weapon 


as x-rays before it begins to blow apart under hydrodynamic action 


depends largely on its yield-to-mass ratio and to some extent on other 


construction details. This fraction may range from almost nothing at 


all (or a very small percent) to a significantly more than 80% of the total 
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energy generated (Glasstone, 1962; Brode, 1964 b). 


Before the air has had a chance to re-radiate any of the energy 


deposited by the x-rays, the bulk of this energy is concentrated in a 


relatively small sphere and at a temperature which is typically of the 


order of several million degrees K. There is, however, a small fraction 


of x-rays, from the high frequency end of the spectral distribution function, 


which penetrates to a distance of perhaps a few meters, and heats this 


shell to temperatures in the 10 I OOOoK range. Energywise this heating is 


insignificant but it makes a contribution to the fireball phenomenology 


which is of some interest. By consulting the table of mean free paths on 
* 
p. 447 of (3) 
we discover that this shell is opaque. As long a s it exists 


such an opaque shell hides the much hotter sphere on its inside and all 


that can be observed is the radiation from the shell itself, which is 


comparatively dim. 


This phase is always very short-lived and terminates when the 


radiation from the center floods into the shell and heats it up. During the 


next pha se the fireball can be characterized rather well as an extremely 


high temperature sphere of air surrounding the nuclear source and showing 


a fairly sharp temperature drop at its edge. The interior of this high 


temperature sphere may be at a fairly uniform temperature I and the whole 


may contain quite a large fraction of the nuclear explosion energy in the 


form of heat. Some small fraction always remains in the dense bomb 


vapors, but most of the early pha ses of the fireball development are quite 


independent of the details of the weapon deSign. The subsequent explosion 


and radiation behavior can be derived almost entirely from the properties of 


this hot air. Such a model will be less true in high altitude or space 


* 
DASA-1917-3. 
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environments where the immediate external surroundings fail to contain as 


thoroughly the explosion energy because they lack sufficient opacity or 


optical thickness. 


Throughout the explosion development, radiant energy is emitted by 


the fireball. That fraction which is transmitted by the cold exterior atmosphere 


is called thermal radiation. The rate of energy emission, or radiated power, 


has shape as shown in Fig. 2-1. If the opaque fringe layer has been penetrated 


early enough and if the instruments used for measuring the thermal power 


have sufficient time resolutions, the signal will include the early peak which 


is shown a's a dotted line. Otherwise, one sees the two-peak curve which 


is drawn as a full line. The explanation of this curve is an important 


objective of any theory. 


Although it is a rather simple exercise, it is instructive to note the 


rather small size of air volumes required to contain the large amounts of 


energy at the high temperatures created by the absorption of the initial 


flux of x-rays. The following table indicates the radii of spheres for a 


few examples of energy content and temperatures. These temperatures, 


of course, are too high for the air to remain that hot for very long, but in 


the immediate first fractions of a microsecond these radii are representative 


of the sizes and temperatures of the earliest (x-ray) fireballs. 


Size of Spheres of Sea Level Air Necessary 
to Contain 1 KT, 1 MT or 100 MT of Energy at 
Various Uniform Temperatures 


i%i~~~~~~rfeoK 
1 KT 
1 MT 
7 1/2 
3/4 m 
7.5 m 


6 
1 
10 


5 
1 1/4 
12 


4 
1.6 
16 


3 
2.1 
21 
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100 MT 
35 m 


46 


57 


74 


100 


In a very few microseconds I these fireballs would have grown much larger 


and much less hot by the continued diffusion of radiation into the external 


cold air. 


For most considerations these earliest phases of x-ray deposition and 


re-radiation remain both obscure and of little probable importance. When 


the flux of source radiations has been sufficiently intense as to completely 


strip the electrons from the air ions, then that volume of plasma can offer 


only Compton scattering as further resistance to the x-ray flux or as opacity 


to its own re-radiation. The most appropriate physical model for the con- 


tinued expansion of this low emissivity, high energy density region is 


neither by hydrodynamics (which requires relatively long times to accelerate 


masses of gas) nor by radiation diffusion.which presumes many interactions 


over any appreciable temperature gradient. The growth of such a heated 


volume is a radiative process which can be characterized roughly by its 


emissivity I temperature I and volume together with the heat capacity of 


the external cold air. The single further physical characteristic necessary 


to include in a growth rate prediction is the fact that the surrounding air 


is essentially opaque to the radiations from this hot air. Detailed 


knowledge of the opacity between this blackness at cold temperatures and 


its transparent nature at sufficiently high temperature is at this pOint 


unnecessary. Thus I the rate of energy lost, expressed as a grey-body loss 


rate, is the rate at which energy is deposited in the cold air at the surface 


of the high temperature sphere, viz. I 


dW = 
dt 


41 


2 
4 
4nR aT e 
(2.2-1) 


in which 
e 
is the emis sivity I 
T 
the temperature of the hot isothermal 


sphere I 
R its radius I and 
dW/dt 
the rate of energy change. When 


an appropriate specific heat is introduced, a differential prescription for 


the volume growth and temperature drop results. 


Following this approximation, one can express the rate of growth I 


dR/dt, in the same terms as 


dR = 
dt 


4 
eaT 


Ep 


in which 
E represents the internal energy per unit mass, and 
p 
the 


density of the air just behind th e front, while 
T 
is the inner temperature. 


The usefulness of this approximation in estimating the rate of growth of a 


partially transparent fireball is largely dependent on the accuracy with 


which average or "effective'" interior temperatures I specific energies, 


and emissiv1ties can be chosen. During the most rapid expansion, the 


interior is likely to be considerably non-isothermal, i.e. I the interior 


may be more than twice as hot as the region just behind the front. The 


dependence on the fourth power of the temperature makes this rate quite 


sensitive to such differences. The most uncertain quantity is likely to 


be the effective emissivity, since it represents some average over the 


emitting region, and may also disguise some geometric dependence - not 


all the radiation being emitted radially. Appropriate choices of effective 


emissivity and temperature may make this simple formula appropriate 


for predicting the growth rate during the subsequent radiation diffusion 


phase. 
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(2.2-2) 


The temperature profiles illustrated in Fig. 2 -2 are typical of this 


early radiative growth for a one kiloton sea level burst. The curves 


represent the air temperatures as a function of radius for six selected 


instants in time. The dashed curve indicates the shock temperatures. 


It is the lowest temperature within the fireball at each of these times. 


After about 15 microseconds, the radiation diffusion growth becomes so 


slow that a shock wave begins to form, to compress the newly engulfed 


air and heat it to a temperature substantially below that of the radiatively 


heated inner sphere. With either the early radiative expansion or the 


subsequent adiabatic expansion behind the forming shock front, the inner 


temperatures drop with time in an approximately exponential fashion. 


During this early growth, the power radiated or the thermal radiation to 


pOints outside the fireball is not a significant fraction of the energy it 


contains. 'The time is short, the size is small, its opacities are high, and the 


fireball exterior so well shields the hotter core that the radiation out is 


les s than half a per cent of the available energy .. 


Of course I the radiative properties are influenced by the air density 


as well as by the temperature, and the gradual formation of the shock 


causes an appreciable increase in the air density at the fireball surface 


(as much as tenfold increase at sea level). In the process, the outer 


surface of the fireball passes from a rather diffuse radiation-driven front 


to a sharp, dense shock front. Fig. 2-3 shows some typical early density 


profiles, in which the shock is seen to grow and the fireball interior is 


seen to expand to much less than the external ambient density. 


Reference to the opacities for air as given in Volume 3 will 
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confirm that the shock front at these early densities and temperatures 


is quite opaque. For instance, at the 250 J.Lsec time of Figs. 2-2 and 2-3, 


the emission mean free path for a shock temperature of 30 ,000oK and a 


density of 8 times normal is about 0.01 cm. The fireball will expand to 


much lower temperatures and much larger size before anything behind 


the shock front will become visible. It is during this period that the 


thermal radiation rate decreases toward a minimum and the fireball appears 


to grow dimmer. (Fig. 2-1, before one millis econd.) 


If the fireball growth rate defined in Eq. (2.2- 2) 
is compu ted for 


the earliest time illustrated in the temperature profiles of Fig. 2-2, 


as suming for the moment an emis sivity of unity, the rate is about 


8 
4.4 x 10 
cm/sec. This rate is too high by an order of magnitude in 


comparison with results of the numerical calculation example. The 


calculation showed that the expansion at this 1. 2 microsecond time was 


s till being determined by radiation diffusion. The calculation, however, 


also treated the earliest times by diffusion, and not (as suggested above) 


by transport within a transparent heated region with a radius less than 


one mean free path for the emitted radiation. The appropriate mean free 


path for diffusion is the so-called Rosseland average, hereafter abbreviated 


as Rmfp. The Rmfp is defined as 


in which 
AV 
co 


co 


S. 


.@L 
AV dT dv 
o 
Cd J 
dBv dv 
dT 
o 


is the spectral mean free path I 
the Planck function I 


and J 
dv 
o 
denoting integration over all frequencies. The Rmfp used in 
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(2.2-3) 


the calculation approaches the Compton limit at high temperatures, however, 


and allows the rate of growth to be equally fast. In fact, without special 


consideration for relativistic effects, the diffusion growth can exceed 


even the speed of light. 


At the earliest time illustrated in Figs. 2-2 and 2-3, the fireball 


has grown to more than one mean free path in radius which reduces the 


effectiveness of the inner temperature in driving the continued expansion. 


A more heuristic interpretation of the growth rate formula allows the 


emissivity to be interpreted as a resistance parameter which reduces the 


growth rate to less than the blackbody rate for that central temperature. 


An alternative interpretation treats this efficiency factor as one which 


compensates for the use of the shielded innermost temperature when the 


effective temperature is at some radial position further out and is lower in 


value, i.e., e ~ (To/Ti)4 
where 
To 
is the effective outer temperature 


and 
Ti 
is the screened central temperature. A crude measure of this 


correction and of an appropriate value for this viscosity constant might 


be the ratio of the Rmfp to the radius of the front, i. e ., the reciprocal 


of the number of mean free paths between the radiating interior and the 


front. For the diffusion approximation, such a correction might better 


be expressed in terms of the local temperature gradient as well. 
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The inner temperature of our example calculation at 1.2 micro- 


6 0 
seconds is around 10 
K (Fig. 2-2) and the density 1s still normal 


-3 mI 
3 
(1.29 x 10 
g 
cm) (see Fig. 2-3). The Rmfp is a bit less than one 


meter, while the radius is about 3.2 meters. Taking 
e 
to be 1/3.2 


8 
brings the growth rate down to about 1. 3 x 10 , which is still high 


compared to that for the numerical calculation. The mean free path 


6 0 
decreases rapidly as the temperature falls below 10 
K, however, and 


since the front at 1. 2 ~sec is at around half the interior temperature, 


a more appropriate mean free path might be between 0.92 (the value at 


6 0 
5 0 
10 
K) and 0.12 (the value at 5 x 10 
K). Taking the average of their 


reciprocals, i.e., averaging the opacities, gives about 0.2, so that 


the correction factor, e, becomes 0.2/3.2, and the corrected rate 


7 
becomes'" 3 x 10 
which agrees well with the growth rate at that time 


from the detailed diffu sion calculation. 


The most appropriate specific energy and density values for use 


in the growth rate approximation are those just behind the front of the 


wave, since it is to those conditions that the cold air is to be heated, 


1. e. I it is that heat capacity that will absorb the subsequent radiation 


energy flux. Fig. 2-4 displays the specific energy profiles for this 


one kiloton example for the same time as those of Figs. 2-2 and 2-3. 


It is interesting to test the simple growth rate formula (Eq. 2.2-2) against 


the fireball growth speed that results from the numerical calculations. 


The calculation should show a rate faster than hydrodynamic shock 


growth until the radiation growth has fallen below the speed of hydrodynamic 


motions, and this simple form should show a comparable rate 


until that time I then a much slower rate as the shock wave takes over. 
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Fig. 2-5 compares these rates for the same time period as covered by 


the profiles of Figs. 2 -2, 2-3, and 2-4 and beyond. In these 


comparisons, several approximations are represented by dashed curves, 


while the numerical calculation rate is shown as a solid curve. The rate 


calculated as blackbody at the inner temperature, shown as circled pOints I 


is clearly too high at all times. Even when the lower temperatures of 


the outer edge of the hot region are used to determine a blackbody rate 


(the triangles A of Fig. 2-5) I the rate is at all times toa high. 


When the radiative resistance parameter is represented as the ratio ofthe Rmfp 


to the hot region radius, using the Rmfp evaluated at the hot interior 


temperature I the modified rate is still high at the early times when 


diffusion is still dominant (the square pOints of Fig. 2-5). It drops 


precipitously as the interior cools and becomes opaque at just the times 


when a shock begins to form (at about 10 microseconds in this example). 


Although this approximation is not correct in value, the sharpness of 


the decrease as hydrodynamics takes over can make it a useful indicator 


of the transition onset, and so a reasonable prediction tool. 


The more accurate estimate of the early diffusion growth rate, 


involving the averaged opacity between interior and front, is also more 


subject to error due to the difficulty in judging appropriate front conditions. 


These estimates are indicated in Fig. 2-5 by diamonds. These values 


are closest to the numerical calculations rate of growth at the earliest 


times when diffusion is the dominant mechanism. The earliest profile 


front temperatures are difficult to define because the front is not sharp. 


The rate of growth estimated at these approximate front temperatures with 
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a corresponding resistance parameter leads to the estimates indicated by 


the triangles t> 
in Fig. 2-5. Again, the shock formation times is denoted 


by a sharp drop in the rate estimated in this manner. Both of the blackbody 


rate curves (upper curves of Fig. 2-5) show a fairly sharp drop at shock 


formation time. Such a simple but uncertain formula may be preferable to 


the use of the radiation resistance notion in determining shock formation 


radius and time. Since in this range of temperature and densities, the Rmfp 


decreases with decreasing temperature as about the fourth power of the 


temperature, using the Rmfp as a correction factor then means that the 


adjustment parameter is as sensitive to temperature changes at the black- 


body rate itself. Such critical opacity dependences may provide some sharp 


distinctions in estimates but at the same time present some hazards in 


chOOSing effective temperatures too casually. 


After shock formation, the rate of growth of the fireball should 


follow the shock growth itself until the shock cools to transparency. 


The shock speed for a strong shock is approximately given by 


. 
R s 


(y +l)P 
s 
s 
2p o 


where 
'V = {P /p E )+1 
P 
is the ambient air density and 
r S 
s/j S 
S 
' 
0 


and 
Es 
are shock front values of pressure, density and internal specific 


energy. This approximation is shown in Fig. 2-5 by the symbol <I . For 


the earliest times, the expansion is faster than this shock rate, but at 


later times it corresponds well. 
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(2.2-4) 


Using the particle velocities (u) at the front and the density at 
s 
the front (through the mass conservation relation) provides the relation 


ugP s 
= 


The rate derived from corresponding values of 
and 
for the 


numerical calculation is indicated in Fig. 2-5 by triangles pointing down 


(\7). After nuclear shock catch-up this curve coincides with the solid 


curve for the directly computed rate. 


In the temperature region of interest, a shock can be pictured as a 


(2.2-5) 


sharp gasdynamic jump imbedded in a region of radiation-induced temperature 


vClriation (Fig. 2 -6). The internal structure of this type of wave ha s been 


investigated extensively by Zeldovich (1957), Raizer (1957), and Heaslett 


and Baldwin (1963), to name a few, all of whom employed the equations of 


steady continuum gas dynamics with gray radiative transport. 


The important feature of this picture is the temperature precursor which 


runs ahead of the sharp front. This precursor is created. by the radiation 


from the high temperature region behind the sharp front. One can estimate 


the temperature of the precursor by equating the power radiated by this 


front with the rate of heating in the precursor. In the resulting relation 


= a T 4 
s 


p , u 
and 
ep 
stand for the ambient air density, the shock velocity and 


the internal energy of the air in the precursor. From the latter quantity and 


the equation of state, one can then obtain the temperature of the precursor. 


Using the Hugoniot relations (Section 5.1 of (4)) and a simple analytic fit 
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(2.2-6) 


to the equation of state ,one obtains a t sea level the relation 


(2 .2- 7) 


5 0 
For a shock temperature of 10 
K we calculate a precursor temperature of 


23 I OOOoK and note that a millimeter layer of air at that temperature is opaque. 
5 0 
Up to the time when the shock temperature drops to 10K all the thermal 


radiation comes therefore from the precursor. To make a quantitative 


evaluation of the power radiated during this phase requires a more detailed 


analysis of the radiative transfer problem. Qualitatively one can see that 


the power must decrease with time and this is the decrease following the 


early peak in Fig. (2-1). 


5 0 
As the shock temperature drops below 10K 
I the precursor cool s to 


where it gradually becomes transparent so that the radiation from the shock 


front begins to shine through. 
When this happens the power-time curve goes 


through the minimum which is shown in Fig. 2-1 as the shock precursor 


minimum (SPM). While the shock front gets more and more exposed I the 


power rise because of the exposure is eventually compensated by the 


temperature drop of the shock itself and at that time the power level reaches 


the maximum which is shown on Fig. 2-1 as the shock exposure maximum (SEM). 


During the pha se folloWing this maximum the rate of thermal radiation 


loss from the fireball can be characterized as that from a blackbody sphere 


at the shock front temperature and of radius equal to that of the shock 


radius. Although such a rate describes the fireball emission I the power 


observed at any distance will contain only that fraction which the cold air 
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outside the fireball is capable of transmitting. To a good approximation, 


that fraction can be calculated by assuming a simple cut-off in the trans- 


mitted spectrum. Values of this fraction 
f(T s' 'J c) 
where 
'J c 
is the 


frequency corresponding to a cut-off at 1860A (representing the edge of the 


02 absorption) are shown as functions of the temperature 
(Ts) 
in 


Fig. 2-7). The fraction is evaluated from a tabulation of the Planck radiation 


function and its partial integral by Gilmore (1956). The fraction is defined as 


X c 
3 
X dx 


X 
e -1 
(2.2-8) 


where 
X = h'J /kT 
and hand 
k 
are Planck I sand Boltzmann IS 
c 
C 
s 
constants, respectively (h~ 6.625 x 10-27 erg sec, k~ 1.380 x 10-16 erg/oK). 


During this phase which lasts until the shock temperature has dropped 


, 


to so Iowa value as to make the shock front transparent, the following 


simple expression characterizes the thermal radiation rate for an air burst 


nuclear explosion: 


P R1 471 R2 a T4 f(T 
,'J ) 
(2.2-9) 
s 
s 
s 
c 


in which 
R 
represents the shock radius, 
T 
the temperature, 
a 
s 
s 
the Stefan-Boltzmann constant (5.672 x 10-5 erg/cm2/deg 4/sec) and 


f(Ts''J c) 
the fraction passed by the cold air. 


Unit optical depth for most frequencies grows longer as the shock 


front cools, so that emission from hotter air behind the front begins to 


shine through. The shock front itself becomes fainter and appears to pull 


ahead of the luminous fireball, a phenomenon which is referred to a s the 
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"breakaway" (Gladstone, 1962, Section 2.110). Because the shock has 


been carrying the shock-heated air outwards with its expansion, a rather 


steep gradient in temperature is maintained just behind the front, so that 


a slight increase in unit optical depth exposes higher temperatures but at no 


appreciable decrea se in radius of effective radiating surface. At this time 


the power curve goes through the principal minimum (PMIN) in Fig. 2-1. 


Fig. 2-8 indicates the geometry of fireball temperatures (in cross- 


section) at a time somewhat beyond the time of minimum thermal power. 


While the thermal radiation increases, and while progressively deeper parts 


of the fireball are exposed I the hydrodynamic expansion dominates so that 


the visible or apparent fireball size continues to grow. Eventually, the 


luminous fireball stops expanding and the power output reaches the final 


maximum (FMAX) . 


Throughout this radiative and then hydrodynamic expansion of the 


fireball, right up to the time of minimum light intensity, something less than 


half of one percent of the total yield has been radiated out of the fireball. 


Both integrals of the measured power-time data from tests and of the simple 


expression given above for radiation from the fireball (as determined by shock 


front conditions) lead to an answer close to 0.44%. In the latter integral, 


the properties of the shock front are sufficiently well defined by almost any 


calculation - even those not accounting for radiation transport in the early 


phases, but necessarily taking account of the real gas properties of air. 


(e.g., Brode, 1956a,b). 


Since the air just behind the shock is much hotter and much 


less dense than the air at the front itself (see Figs. 2-2, 2-.3, and 


2-4), the rate of thermal radiation increases rapidly when that air is 


exposed, until the hottest temperatures at the back of the steep gradient 
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behind the shock front become visible and are radiating directly to the 


exterior. Thereafter, as the size of the radiating sphere shrinks and 


the interior cools, the rate decreases. This is the period in which the 


fireball history comes closest to the cooling wave notion expressed in 


a simple form by Zel'dovich, Kompaneets and Raizer (1958) and applied 


into a fireball theory by Bethe (1964). The notion is that a recognizable 


and fixed form cooling wave erodes the hot fireball interior, beginning 


at the exterior and working inwards. After the shock front has become 


transparent, such a cooling wave process is very likely operating, but 


it 1s not at first working into a .f1xed or uniform temperature or density I 


and it is not shrinking the fireball. The outward hydrodynamic expansion 


is still too strong. When the outer regions have all become sufficiently 


cool and transparent so that the inner radiation-heated region is exposed, 


then the conditions suggested for a cooling wave are approximated. Even 


then the temperatures are not constant and the surface area is shrinking 


rapidly I so that the cooling rate decreases. When this interior sphere 


has cooled to below about 10 I OOOoK, the whole of the fireball has become 


relatively transparent, and the subsequent radiation losses are characterized 


more by a grey body approximation, i.e. I characteristic of a volume of air 


of low emissivity - one of less than unit optical thickness .. It may also 


still be expanding adiabatically I and contributing energy to the shock 


growth. 


Temperature profiles spanning this period from principal minimum through 


final maximum and on to a transparent fireball are illustrated in Fig. 2-9. 


For a yield of one kiloton I the cooling wave is les s obvious as a wave 


than as a rather sudden depletion of the hottest interior region. At larger 


53 


yields, where more optical thickness is represented at every stage, 


the progress of a cooling wave from outside toward the center is more 


easily imagined (Brode, 1964a, Fig. 5, b Fig. 15). 


In this rather complex power radiated history of two or three maxima, 


as illustrated in Fig. 2-1, the final pulse represents a total energy of 30 


or 40% of the total yield of the nuclear device. When all the energy is 


accounted for, including that in the infrared which originates in shock 


heated air outside the visible fireball and is radiated only very slowly, the 


fraction may be even larger. 


There are several features of this one kiloton explosion that have 


not yet been mentioned and that are of lesser influence on the thermal 


radiation and fireball behavior at sea level, but which become relatively 


more important at other yields or altitudes. One such feature is a second 


shock wave which originates within the bomb vapors, traverses the early 


sphere of hot air behind the radiation front, and overtakes the strong 


shock that forms the fireball surface at later times. This debris or 


bomb shock is seldom in evidence in sea level explosions, and has lost 


most of its energy long before it overtakes the main shock, so that it 


contributes little to the fireball surface or thermal radiation histories. 


Because the hot interior of the fireball is for most of the fireball expansion 


a region of long mean free path, it is a region of nearly uniform temperature. 


When the case shock compresses and heats this air further, some of that 


heat is promptly re-radiated ahead, forcing this interior shock to behave 


isothermally rather than adiabatically. This isothermal shock can lose 


energy very rapidly by this means, and may persist only through the 
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continuation of its outward momentum. 


A history of the radii of this shock and other fronts in this kiloton 


example is shown in Fig.2-10. When this debris shock travels outward 


to the edge of the fireball, it encounters a sharp discontinuity in density. 


At .that pOint, a reflected shock originates and is returned inward to 


implode upon the Origin. Here again, is a phenomenon which has no 


consequence for this example, but may be prominent in high altitude events. 


The vaporized bomb expands along behind this debris shock,-' but at sea 


level is not visible until very late - after the second maximum in the 


thermal power. This bomb debris is not realistically treated in any of 


the usual calculations I since they invariably assume radial symmetry 


and allow no mixing or turbulent flow. When it emerges in the transparent 


fireball at late times, the vaporous debris has become highly turbulent 


and has evidently mixed with considerable fireball air. 


Although Fig. 2-10indicates a transition from radiation expansion 


to strong shock expansion, the radiation diffusion does not stop. As the 


shock brings down the density in the interior air, the opacity of that air 


decreases also, and the radiation is allowed to diffuse into some of the 


now shock-heated air. The dotted curve below the shock front curve of 


Fig. 2-1 0 indicates the position of the radiation front. Most of its outward 


excursion is due to the flow of air in the expansion behind the shock 


itself. At times later than shown in Fig. 2-10, the radiation front and the 


visible fireball drop behind. The short dashed curve near the end of the 


shock front curve of Fig. 2-10represents a position close to the fireball 


front - being the locus of points at 50000 K - with higher temperatures 


inside of that radius, and colder temperatures outside. 
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The continued flow of radiation is made more obvious in a plot of 


the temperature histories of several shells of gas representing the air 


that was shocked to a particular temperature, cooled in the subsequent 


adiabatic expansion, but then reheated by the radiation wave following. 


Such a set of curves are shay rL in Fig. 2-11, where at particles shocked 
5 
a 
to 10 , 70,000 and 40,000 K the adiabatic cooling is arrested by the 


arrival of the radiation diffusion wave which causes that shell of air 


o 
to rise in temperature again. The air starting at the 20, 000 K shock 


pOint is never over-run by the radiation wave, i. e., the radiation wave 


stops before it gets that far, having run out of energy and not being 


aided by further expansion which would help to reduce the opacity of the 


cooler air in front of it. 


A great many nuclear weapon applications, tests, and effects 


interests involve the thermal and fireball effects of nuclear explosions on 


or close to the surface of the earth. -Many interesting and novel inter- 


actions occur which are not evident in air bursts well away from the 


sur face. However, there is no intention of providing a review of these 


factors in this current effort. It should suffice to point out that all of 


the essential features which are described and followed here are also 


an important part of surface bursts, while the latter are further complicated 


by the early injection into the fireball of massive amounts of earth material, 


and by the geometric distortions of the fireball that occur as a consequence 


of shock and thermal reflections from the earth I s surface. The change in 


radiator shape from spherical to at best hemispherical or worse a partially 


obscured hemisphere means that the thermal flux to other points on the 


earth1s surface will be less than that from an air burst. Total flux at 
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pOints in the air above the burst may at the same time be increased. 


As the earliest pictures of nuclear explosions (G1asstone, 1962) 


clearly show a further consequence of the ground involvement is the 


"dust skirt- which precedes the fireball shock and largely obscures 


the base of the fireball. Although not visible in any of the pictures, 


there must also be vast amounts of earth shovelled into the hot fireball 


interior at an early time (Brode and Bjork, 1960), and this rna teria1 cannot 


fail to have profound effects on both the temperature and thermodynamic 


state of the fireball gases and on the opacities or optical properties of 


that region. Test observations .indirectly attest to the influence of such 


surface effects. 


Observations and measurements at very late times in the fireball 


history show that the radiation rate trails off with a very long tail (as in 


Fig. 2-1) and comes from shapes other than simple spheres. The fireball 


at late times is like a bubble in the atmosphere - having very low densities 


in its interior - and so it rises, and in rising breaks up at the bottom to 


transform itself into the familiar atomic cloud ring or torroid which rolls 


its way up through the atmosphere. The torroidal circulation that is 


induced is quite strong and serves to severely limit mixing of the hot 


fireball gases with the exterior cold air, thus prolonging the existence 


of air and debris at temperatures of thousands of degrees Kelvin, while 


the cloud rises in the atmosphere. When much earth material and/or 


water vapor is present, the late fireball remains opaque, and the rate of 


late radiation is more determined by the rate of turbulent mixing which 


brings hot gases to the cloud surface rather than by the radiation transport 
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properties alone. For an air burs t well above the surface, however, the 


late fireball becomes quite transparent, so that only a faintly luminous 


ring assures us that the rise and circulations are much the same as for 


lower burs ts • 


2.3 
Other yields and altitudes 


The example of a one kiloton detonation at sea level contains all the 


basic physical phenomena which enter into consideration at [;ther yields 


and at altitudes up to about 70 km. The overall appearance is I nevertheless I 


appreciably different' since the individual events which are responsible for 


the various maximum and minima in Fig. 2-1 occur at different times. 


In carrying out a discussion of these changes I it is useful to note 


that the relation between shock radius and time can be approximately repre- 


sented by a hydrodynamic scaling law. To formulate this we introduce the 


scaled variables 


where 
Y 
is the yield of the explosion and p = pip 0 
the a mbient air 


density relative to that at sea level. In our 1 kiloton sea level exa mple the 


scaling factors are of course unity. 
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(2 .3 -1) 


(2.3-2) 


The scaling law is not valid until after the debris shock has caught 


up with the somewhat slower shock which is driven by only a fraction of the 


total yield. The scaling law can be deduced from the strong- shock solution 


for a point source (Taylor, 1950; and Sedov, 1959). This limits the validity 


of. the scaling laws at late times when the shock becomes weak. The law 


obtained in this manner takes the form 


* 
* 2/5 
Rs == k(t ) 
(2.3-3) 


where the subscript 
s 
denotes that the value is taken at the shock front. 


From the radius 
R 
== 20 m 
read off Fig. 2-10 for 
t = 1 msec, 
s 


one determines the proportionality factor to be 


k = 20 m (msec)-2/5 
(KT)-1/5 
(2 .3-4) 


In the above scaling law the scaling factors cancel out of the expression 


for the shock velocity 


dR s 
= 
* 
dt 


* 
= 2 
5 


* 


* -3/5 
k(t ) 
(2.3-5) 


which makes this velocity a function of 
t 
only. Applying the Hugoniot 


relations one can now show that the temperature 
Ts 
behind the shock is 


S9 


* 
also very nearly a function of 
t 
only. This is not an exact result because 


it depends on certain assumptions about the equation of state which are only 


approximately true. If one checks the prediction that 
T sis a function of 


the scaled time only against computed results, one finds that it fits the 


changes with yield at a given altitude very well. The changes with altitude 


at a given yield are not gi ven with quite the same accuracy, but are still 


suffiCiently close for most purposes. 


It should be noted that the above scaling procedure differs somewhat 


from the so-called Sachs scaling where one introduces the variables 


(2.3-6) 


(2.3-7) 


If one expresses the ambient pressure and density 
p 
and 
p , the yield 


Y 
and the variables 
Rand 
t 
all in the same system of units/the scaled 


variables are dimensionless. It is more convenient, however, to replace 
p 


and 
p 
by the ratios 
15 = plpo 
and is = pip 0 
relative to the sea level 


values, and to express 
Y as before in KT. With this choice, the strong 


shock relation between Rand t 
is the same as between 
R 


i. e. Eq. (2.3-3) with the same value of the constant 
k . 


* 
* 
and 
t 


The two methods of scaling differ in regard to what are considered 


similar situations. For the starred variables similarity implies that for example 


the hydrodynamic velocity and the temperature are unchanged; for the variables 


with the tilde the Mach number and the temperature ratio TIT 
- 
0 
are unchanged. 


Either choice is acceptable, but ours has the advantage of using only one 


parameter to characterize the altitude. 
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From the time of shock formation until breakaway the thermal radiation 


comes partly from the shock precursor and partly from the shock front, and 


it is evident that the shock temperature is a major factor in determining the 


timing of the maxima and minima during this period. At a given altitude where 


one ha s a one-to-one relation between shock and precursor temperature 


(see Eq. 2.2-7 for the sea level case) it is fairly accurate to state that the 


shock formation maximum, the shock precursor minimum and the shock exposure 


maximum occur at fixed values of the shock temperature and therefore at 


fixed values of the scaled time. As one considers different altitudes the 


relation between 
Tp 
and 
Ts 
changes and one finds different values 


of the scaled times associated with these features of the power curve. 


After breaka way the radiation comes from points to the inside of the 


shock front whose locations depend on the optical properties of the air and in 


turn on th~ temperature and density distribution. This is a radiative transfer 


problem and hydrodynamic scaling, where times vary as the cube root of 
Y , 


is replaced by radiative scaling, where times vary approximately as the square 


root of 
Y 
(Glasstone, 1962, section 7.92) . 


Altitude scaling is a more difficult problem than yield scaling. We have 


already mentioned the effect of the changing relation between 
T sand Tp 


To this we must add that the relative importance of hydrodynamics and 


radiation transfer shifts with increasing altitude in favor of the latter. Thus 


shocks form more slowly and radiation is emitted more rapidly a s one goes to 


higher altitudes. As a result the features before breakaway are increa singly 


delayed and the maxima and minima tend to become weaker. The final 


radiative pulse on the other hand advances in time and becomes more prominent. 


At about 50 km the early features have become washed out and what was the 


final pulse is now the only pulse. 
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Chapter 3. 
SUMMARY 


As the reader of thi s report will have gathered any attempt at following 


the evolution of a fireball by numerical means utilizes a whole spectrum of 


facts and assumptions ranging all the way from being undisputable to being 


highly suspect. This is likely to leave him with a somewhat uncomfortable 


feeling about the reliability of such a calculation. In this summary we will 


put the finger on some of the underlying assumptions I point out what we 


know about their validity and evaluate how strongly our lack of basic 


information or of the willingness to spend computing dollars will influence 


the final product. 


3.1 
Equation of state 


Nearly all calculations make the basic assumption that the air remains 
* 
throughout in a state of LTE. Once this is accepted it follows that the 


relation between the various state variables can be found by the methods 


of statistical mechanics. The application of these methods is very straight- 
** 
forward and the results as presented in (1) and (3) 
are probably correct to 


within a few percent. In some instances analytic fits which were made to 


feed these results into a computer have been poor but this problem can 


certainly be overcome and should not contribu te significantly to errors in 


hydrodynamics or other phases of the main calculation. Some problems may 


arise in the central region where one has debris rather than air and even 


more so in the transition region where one may have a debris air mixture. 


Fortunately many important results are rather insensitive to these details. 


* Local Thermodynamic Equilibrium (LTE) . 


** (1) stands for DASA-1917-1, etc. 
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3.2 
Absorption coefficients 


This subject has been discus sed in detail in (2) and related facts 


are brought up in Chapter 1 of this volume. There are several ways of 


describing the absorption which differ in the amount of detail which is 


presented. The most detailed description consists of a listing of lines 


with intensities and line shape parameters on top of a continuum. All 


these factors are subject to errors as we shall briefly discuss. 


In the low temperature case where the lines are due to molecular 


systems the information comes largely from experimental spectroscopic 


studies. The limitations of our knowledge about frequencies and intensities 


is discussed in (2) Chapter 7. The information on line shapes is almost 


non-existent. To this one should add that one can hardly afford to include 


any but the strongest band systems. Even the rather minimal choice of 


eight band 'systems in the most recent version of the SACHA program brings 


the number of transitions to over 190 I 000. 


In the high temperature case the absorption comes from inverse 


Bremsstrahlung and from transitions in atoms and atomic ions. There is 


a strong continuum due to the first and due to photoionization. There is 


a fairly well developed theory and some experimental information backing 


it up. On top of the continuum is a large number of lines. A few levels 


have been observed experimentally but the majority I especially for the 


highly ionized atoms I have not been observed and must be obtained 


theoretically. . It is certainly neces sary to find the transition probabilities 


by quantum mechanical methods. These are so complex that one is forced 


to make radical approximations to get any answers and the results are not very 


reliable. 
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The calculation of the line contribution is the most elaborate part 


of the program and again one can hardly afford to include any but the 


strongest lines. This involves a somewhat arbitrary cut-off procedure 


whose practical effect can only be evaluated when one specifies how 


the absorption coefficient is to be used. 


The detailed description of absorption with fine spectral resolution 


greatly exceeds the requirements of radiative transfer calculations. As 


shown in (2) Chapter 2, it is unfortunately difficult to define averages which 


permit satisfactory calculations. Thus Planck and Rosseland means which 


and 
-1 
average 
~~ 
~~ 
respectively apply only in limiting situations; 


the one for very transparent, the other for very opaque media. Nevertheles s 


such means are useful and have been calculated. In the specific case of 


line effects, mentioned in the preceding paragraph, the contribution is 


not very large until one reaches temperatures like 2 x 10 5 OK 
and high 


densities. 


Because of the many uncertainties entering the calculation of 


absorption coefficients one has no systematic way of estimating their 


accuracy. The responsible authors of opacity calculations are generally 


confident that their results lie within a factor of three of the true values. 


In the intermediate temperature range where the opacity reaches a 


maximum ,the accuracy is probably somewhat better. Because of the large 


opacity the averaging procedures I which are appropriate for radiative 


transfer calculations, put the most weight on those parts of the spectrum 


where 
is small and very little weight on the lines. The Rosseland 


mean does and the Planck mean does not fall into this class. Because of 
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the emphasis on the continuum, where one has more reliable information, 


the Rosseland mean is expected to be more accurate. 


The extent to which opacity errors falsify fireball calculations 


depends on the temperature range. Inspecting the temperature profiles 


of Fig. 2- 2 which are typical for the early stage of a fireball one finds 


large temperatures near the center which makes the air in that region very 


transparent. Further out the temperature drops so that the opacity rises, 


goes through a maximum and then drops again. In the transparent central 


region radiative heat transfer is rapid and keeps that region at a fairly 


uniform temperature, as one sees in Fig. 2- 2. Just how uniform this 


profile is has very little effect on the rate of expansion and therefore 


opacity errors by a factor twice or even more are not serious in that region. 


The opaque zone around the central region acts as a radiative barrier 


and the development of the fireball does depend quite critically on the opacity 


there. During the very early phase where hydrodynamic motion is still 


negligible compared to the radiative expansion the section of the opacity- 


temperature relation near the maximum determines the rate of that expansion. 


It also determines when and where the hydrodynamic shock begins to form. 


When shock temperatures are still high, the opaque shell forms in 


a temperature toe ahead of the shock. This is the precursor which has been 


sketched in Fig. 2-7 and which causes the early structure in Fig. 2-1. 


At this stage the opacity is still of interest, since it determines the 


character of the escaping radiation and other observable phenomena, but 


the rate at which the fireball expands is given by the shock speed which 


does not depend on the opacity in the toe. 
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The next phase starts when the shock temperature has dropped low 


enough that the shocked air becomes opaque. This is aided by the high 


density directly behind the shock which can be seen for example, in the 


density profiles of Fig. 2 -3. Up to that time radiative transfer plays a 


major role in feeding energy to the expanding shock front. Now that 


source fades out and hydrodynamics takes over as the dominant mechanism 


for energy transfer. The details of the. change depend quite critically 


on the opacity relation. 


Upon further cooling the shock front becomes transparent again 


and the opaque shell recedes toward the center. This starts the long time 


interval during which the power vs. time curve of Fig. 2 -1 goes through its 


minimum, rises back to the final 
maximum I and starts to drop again. The 


calculation of this phase also depends quite critically on the opacity. A 


test calculation made with an opacity twice the accepted value stretched 


the total duration of this phase by almost a factor of two with a corresponding 


reduction of the maximum power level to about half of what it was in the 


earlier calculation. Thus, errors in the opacity relation could lead to fairly 


severe discrepancies between fireball models at the time of the second 


maximum. 


3.3 
Radiation hydrodynamics codes 


The purely hydrodynamic part of any code is probably as accurate 


as the equation of state that is being used except for the smearing out of 


the shock front introduced by the artificial viscosity method. The accuracy 


of radiative transfer calculations is less certain, unless one is justified 


in using the diffusion approximation. In that case the limiting factor is 
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probably the accuracy of the opacity. Difficulties do arise I however I at 


the front of an opaque shell. Consider I for example I two zones labeled 


a 
and 
b 
whose temperatures place them on the rising branch of the 


opacity curve as indicated in Fig. 
3- 10 
As the heating wave progresses 


thepe points move up on the curve. The radiation escaping to the outside 


comes at first from zone 
a 
and passes without attenuation through zone 


b 
. As zone 
a 
climbs higher the radiation rises as 
T4 
but when zone 


b 
becomes sufficiently opaque to take over the role of zone 
a 
the power 


output drops. This cycle is repeated when zone 
c 
and others after it 


climb into the position originally held by zone a. The result of this is 


a sequence of maxima and minima in the power versus time curve which 


has no physical reality. This spuriou s effect can be counteracted by using 


finer zone sizes but at the expense of increasing the running time which 


increases as the square of the number of subdivisions per zone. Actually 


this is not necessary I since test calculations show that the cruder zone 


divisions lead to the same average power and to the same rate of expansion 


as a very fine division. A related problem arises when the artificial 


viscosity routine introduces improper heating ahead of the shock front. 


Letting the pOint 
a 
in Fig. 
3 -1 represent the shocked zone and band 


c 
the zones just ahead of the shock I this heating would make points 
b 


and 
c 
lie at too high temperatures. The calculated attenuation of the 


radiation from the shock is therefore larger than what it should really be 


and leads us to predict too low a brightness of the fireball. The reduced 


output has I however I practically no effect on the calculated motion of the 


fireball air because at that stage the amount of energy lost by radiation 


is still too small to influence the hydrodynamics. 
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Other errors may be introduced by the use of approximate integration 


routines such as the multiple ray or the two stream techniques which have 


been discussed in Chapter 1. Given a set of experimental data one can, 


within limits, adjust the opacity temperature relation so that either model 


will reproduce these data. It is therefore not really possible to disentangle 


errors which may arise from the use of these models on the one hand, and 


from incorrect absorption coefficients on the other hand. 


3.4 
Deviations from LTE 


As we have repeatedly stated, nearly all calculations assume the air 


to be always in LTE. There are, however, some equilibration processes which 


are decidedly slow on the time scale of nuclear fireballs. At somewhat 


elevated altitudes one finds for example, that the proces ses responsible 


for populating the vibration ally excited levels of 
02 
and 
N2 fall into 


this class. These processes are discussed under the heading "vibrational 


relaxation" in (4) section 6.1 In the case of 
02 ' populating the 


vibrational levels reduces the photon energy required for reaching the 


Schumann-Runge continuum below the 8.5 eV which it takes from the ground 


state. As long as these levels are not populated the actual absorption is 


therefore less than it would be in equilibrium. Similar considerations apply 


to the Birge-Hopfield transitions in 
NZ . In most codes these delays are 


just ignored. Hillendahl (see Appendix A) has attempted to account for them 


by means of a fairly crude model assumption. 


Other deviations from LTE are caused by the slowness of chemical 


reactions at temperatures I say I below 6000 0 K (see (4) sections 6.9 


and 7.5.) Among the molecules which can form in this temperature range is 
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N02 
which ha s a large absorption coefficient. The delay in forming this 


molecule when air is suddenly heated by a shock and the subsequent delay 


in removing it again when the air cools down can change the absorption 


significantly from t he value at equilibrium conditions. If the temperature 


drop is rapid enough the 
N02 
concentration may stay for a long time at 


the high concentration corresponding to 3000 0 K even though the temperature 


o 
has dropped below 2000 K. In this situation one speaks of 
N02 
as being 


frozen in. 


Non-equilibrium processes also occur at the debris air interface. It 


has been pointed out in (4) section 5.2 that this is very poorly understood. 


It is, in particular, quite uncertain what temperature the shocked air 


would reach and what X-ray spectrum would be emitted by that air. 


The questions raised in this chapter clearly do not exhaust the subject 


of possible errors in the present state of the art. It is, in fact, quite likely 


that effects with more practical significance have been overlooked. Still, 


this enumeration should provide the reader with some guidance what he 


should watch. 
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T 


FIG. 
3-1 LOCATION OF 3 ADJACENT ZONES 
ON OPACITY CURVE 
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Appendix A. 
A Radiation-Hydrodynamics Code 


A . 1 
Introduction 


In this appendix, a sample radiation-hydrodynamics code is 


presented which employs, with yarying degrees of sophistication, much 


of the physics and basic data presented elsewhere in this volume. 


In keeping with Chapter 2, this code describes the radiative and 


hydrodynamic properties of a sphere of hot air. Details of the weapon 


itself are not of interest in the present context, and rather crude 


generalizations have been used to represent the gross properties of the 


hardware. 


The code is presented as a means of demons trating some of the 


techniques of radiation-hydrodynamics, as described in Chapter I, the 


application of basic physical data, as described in Volume 2 , 


and as an Hlustration of the results discussed in Chapter 2. The code 


is not intended as a demonstration of the programming art and has not 


been polished-up for presentation here. A great deal of the program 


could be deleted were the program to be used only for present purposes. 


Much of the basic philosophy of this code has been presented in Chapter 1 


and by Hillendahl (1964), and will not be repeated in detail. 


The basic equations of the problem are the conservation equations 


of radiation-hydrodynamics for a one-dimensional spherical system 


which can be written in Lagrangian form as 


au = 
at 
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Conservation of momentum 
(A. 1) 


u = aR 
Defini tion of velocity 
(A.2) 
at 


v = 4 R2 a.R 
n 
am 
Conservation of mass 
(A.3) 


Q = ~[,2 
on compression 
c 
V 
at 
Definition of Artificial 
(A.4) 
Viscosity 


Q = 0 
otherwise 


Conservation of energy 
(A. 5) 


E = E(V, T) 
(A.6) 


p 
:: 
P(V, T) 
(A. 7) 


:1 
is an integral functional of V 
and 
T 
(A. 8) 


where 


U 
= local fluid velocity 


t 
- 
time 


R 
= radius 


p 
= pressure 


Q 
:= artificial viscosity 


m 
= mass 


c 
z an arbitrary constant near unity 


po· initial density 
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v 
= specific volume (reciprocal density) 


:1 
= radiative net flux at R 


E 
= internal energy 


Quite generally, the mathematical formulation of the problem may be 


characterized as an initial value problem whose solution consists of a 


time-wise and mass-wise integration of a well defined set of hyperbolic 


partial differential and partial integro-differential equations. 


The solution of these equations is carried out by numerical techniqut~s 


in which values of the dependent variables are determined in terms of the 


two independent variables (the time and lagrangian zone mass) by means of 


finite difference equations which are used to represent Eqs. (A. 1) - 
(A. 8) • 


For purposes of numerical computation, the fireball configuration is 


represented by a series of conc?ntric, contiguous, spherical mas s shells. 


The mass of the k th zone is designated by 
mk (gm/ cm3). Since the mas s 


zones retain their identity throughout the time-wise development of the 


configuration, the zone index 
k 
and the time 
t 
(secon ds) are convenient 


choices for the independent variables. 


Integration of the set of 8 equations (Eqs. (A. 1) - 
(A. 8)} then deter- 


mines the values of the 8 dependent variables as functions of 
k 
and t. 
-1 
U(k, t) (cm sec 
) and 
R(k, t) (cm) are used to specify the ins tantaneou s 


values of the interface velocity and radius of the outer surface of the k th 


mass shell. 
:1(k,t) (ergs cm- 2 sec-i) is used to specify the instantaneous 


value of the net radiative flux at the ou ter boundary of the k th mas s zone. 


P(k,t) (dynes cm- 2), Q(k,t) (dynes cm- 2), V(k,t) (cm3 gm- 1), T(k,t) 


(oK), E(k, t) (ergs gm -1) are used to represent the instantaneous values of 
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the pressure I artificial viscous pressure, specific volume (reciprocal 


density), temperature and internal energy of the mass zone 
ffik • 


In a purely hydrodynamic problem without radiative transfer, it is 


the standard practice to reckon the thermodynamic properties of a zone 


(i.e.: pressure, internal energy, density, temperature) as constant 


average values over each zone. These values are also considered to 


be the central values of these variables at the geometri cal zone centers. 


Particle velocities are reckoned at the zone interfaces; the interface 


density and pressure gradient are formulated in terms of the values at 


the zone centers. The zoning mesh must be chosen fine enough so that 


the variation in properties from zone to zone is small enough to insure 


that these average values are meaningful. 


I n a problem which also includes radiative transfer, the above 


restrictions must also be met. Zone sizes in a problem including radiation 


will generally be smaller than the zone sizes required by hydrodynamics 


alone. The addition of radiative transfer to the problem will, in general, 


add further restrictions. 


If the temperature is taken as cons tant acros s each zone, temperature 


discontinuities will occur at the zone interfaces.' Radiative variables like 


T4 
will have even greater discontinuities. More detailed examination 


indicates that the temperature and its spatial derivatives should be 


continuous at the zone interfaces. Thus, consistent with the expansion 


used in Eq. (2.5-9), the source function 
B is taken as linear between 


zone centers. Then the discontinuous spatial derivatives of the source 


function which occur at the zone centers do not appear in the formulation. 


In a more general formulation I a higher order polynomial could be used to 
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fit the source function through the zone centers, but numerical experience 


has indicated such a procedure resulted in only minor improvement in the 


computations. 


The central zone temperature 
T 
is used as the average over the 


zone for purposes of computing average zone pres sures, internal energies, 


and 'is used also in the 
Z, A, and 
W 
computations (Eqs. (2.5-13) 


through (2.5-15». This is done primarily for purposes of convenience, 


but can be at least partially justified. In regions of small temperature 


gradient, no problem occurs since the central and average zone tempera- 


tures are nearly identical. In regions of large temperature gradients I the 


"average value" of the temperature is poorly defined in terms of the rapidly 


varying radiative variables, and it is preferable to keep the problem well 


pOised hydrodynamically. 


The specific volume 
V 
is taken as having a linear variation 


, 


acros s each zone. The values of 
V 
at the zone centers and zone 


interfaces are then uniquely defined and afford no further difficulty. 


The 
Z, A, and 
W 
functions are then computed using the average 


zone temperature 
T 
and specific volume V. Use of the average 


specific volume is justified since these functions show a relatively weak 


density dependence.- 
Neglecting the variations in temperature across 


the zone causes only relatively small errors in the high temperature inner 


fireball regions since these functions show only a weak dependence upon 


temperature. In the low temperature regions I the 
Z, A I and 
W 


functions vary about as the ninth power of the temperature. Even with 


extremely fine zoning I large temperature gradients occur across each 


zone, and the 
Z, A I and 
W 
functions would be ill defined in terms 
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of any average temperature no matter how the average be defined. 
Bu l 


the emission from these low temperature regions is small compared to 


the emission from the high temperature regions further inside, and these 


layers act primarily as a selective absorber for radiation from larger 


optical depth s • Hence the 
T 
Z 
function must be known with some 


accuracy I but small errors in the 
A 
function, originating because of 


the use of the average zone temperature I can be tolerated. 


The 
Z+ 
function I however I has peculiar properties in the low 


temperature region which allow its values to be obtained with sufficient 


accuracy. As discus sed in Chapter 8 I the spectral absorption coefficient 


varies extremely rapidly with wavelength so that the spectrum is effective. y 


divided at some wavelength into absorbed and transmitted fractions. This 


transition wavelength I however I depends only weakly on the zone tempera- 


ture I and hence the average zone temperature will again suffice. 


It should always be born in mind that these numerical approximations 


all improve as the zone size is decreased and can, in principle be made 


accurate to any desired precision. In practice I however, the finenes s of 


the zone mesh is limited by the cost of computation. Skill is thus required 


to accomplish a large computational program within a limited budget. One 


tries to test each situation for sensitivity to zone sizes and achieve a 


compromise between economy of computation and accurate representation of 


the physics. 


For purposes of carrying out the integration procedure by numerical 


methods I the basic equations (Eqs. {A. 1) - 
(A. 8» are replaced by a set 


of centered finite difference equations as follows. The notation and 
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centering can best be seen by reference to Fig. A-I. 


n+l 


°k-l·/2 = ( 
c 
~ ~ n\ 
)(v~~i/2 -V~ _1/2)2 
n+l 
n 
2 
n+l 
n 
Vk - 1/ 2 + Vk - 1/2 
4n R~+l 
t 
- 
t 


n+l 
= 
0k-1/2 
0 
if 


Ek- 1/ 2 - Ek - 1/2 
+ Pk- 1/2 + Pk - 1/2 + On+ 1/2 
k-1/2 - Vk - 1/2 
n-t 1 
n 
~ n+ 1 
n 
~n+ 1 
n j 


tnT 1 _ tn 
2 
k -1/2 
tnT 1 _ tn 


+ 1u.ifRn+1)2 9 n+ 1 _fRn+1\2 9 n+ l 
+ IRn\2 9 n - 
fRn Y 9 n I 
= 0 


mk 
\ k 
k 
\ k - Y k - 1 
\ k J 
k 
\ k - 1 J 
k - 1 
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(A.9) 


(A. 10) 


(A. 11) 


(A. 12) 


(A. 13) 


n+ 1 
( . .n+ 1 
n+ 1 ) 
Ek- 1/ 2 ~ E vk- 1/ 2 , Tk- 1/ 2 


:1 n+ 1 = 9. fvr:+ 1 Tr:+ ~ 
i = I, ••• •••••••• L 


k 
k~ 1 
, 
1 J 


Eq. (A.9) expresses the conservation of momentum and is centered 


over the grid point (n, k) • 


Eq. (A. 10) simply expres ses the definition of velocity and is 
, 


centered ov~r the grid point (n-ti-l/2 ,k). 


The conservation of mass, expressed by Eq. (A. 11), is centered 


over the grid point (nT 1/2 ,k -1/2, i. e. over 
Q~: ~;~ . 


The difference expression (Eq. (A. 12» for the artificial viscosity 


is centered on the grid location (n+ 1/2 ,k-l/2) except for the 
R~+ 1 


factor. The artificial viscosity is thus correctly centered for use in the 


energy Eq. (A. 13), but lags a half time step behind in the equation of 


motion (Eq. (A. 9). 


The energy equation (Eq. (A. 13» is centered over the grid location 


(n+1/2 ,k-1/2) and constitutes an implicit expression to determine the 


local temperature 
n+1 


Tk - 1/ 2 • 


Eqs. (A. 14) and (A. 15) are equations of condition rather than finite 


difference equations, and expres s the equations of s tate for the fluid at 


the location (n+ I, k-l/2) in terms of 
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and 
nT 1 
Tk - 1/ 2 • 


(A. 1.4) 


(A • 15) 


(A. 16) 


Eq. (A. 16) , which is written above only in symbolic form, expresses 


that the local net flux .1~+ 1 
i~i an instantaneous integral functional of the 


temperatures and densities of all of the L zones in the configuration. The 


set of Eqs. (A.9) - 
(A. 16) thus constitutes a set of 8L equations in 8L 


unknowns. 


The following set of 11 auxiliary equations are used to evaluate 


Eq. (A. 16): 


n+l 
{n+I 
n+I 
n+l 
n+J 
ZOk-1/2 = ZO\k~l/2' Tk- l / 2, TRk- l / 2, t 
/ 


n+l 
ZIk - 1/ 2 


n+l 


AHk - l / 2 


= ZI(n+ 1 
Tn+ 1 
tn~ ~ 
k-l/2' k-l/2' 
/ 


( 
n+l 
n-lrl 
n+~ 


= AH Vk - l / 2 ' Tk- l / 2 , t 
) 


n+l 
(noTl \ 
Ak-l/2 = exp - 6.T k-l/1 
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(A. 17) 


(A. 18) 


(A. 19) 


(A. 20) 


(A. 21) 


(A.22) 


Sn+ 1/2 
k-1/2 
= 


W,1+ 1 
= 
k-l/2 


9 n+ 1 


k 


n+ 1 
n+ 1 
BCk- 1!2 - 
BCk+1!2 


n+ 1 
n+ 1 
6T k_1/ 2 + 6T k+ 1/ 2 


n+ 1 
n+ 1 
n+ 1 
1 - Ak-l/2 - Ak-1/2 6T k - 1/ 2 


zepn+ 1 + Fln+ 1 [1 _~R~: i)2] ZIn+ 1 
k 
k 
R~l 
k 
k 


where 
(J = 5.67 x 10- 5 is the Stephan-Boltzmann constant, and 
TR 
is 


the temperature used to specify the spectral character of the radiation 


incident upon the zone from regions of smaller radius. 
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(A. 23) 


(A.24) 


(A.25) 


(A. 26) 


(A. 27) 


(A. 28) 


The analytic expressions for Eqs. (A. 18) and (A. 19) are not given 


explicitly in the written text, but are available in the FIRE BALL code 


listing which follows. 


In carrying out the numerical integration scheme, it is assumed 


that the initial values of 
R, T, V, E, P, U, Q, and :1 are known at 


the initial time 
tn 
for all values of 
k • The difference equations 


(Eqs. (A. 9) - 
(A. 28» are then used to determine values for these variables 


"1 
at 
tn 
= tn + Li t 
for all values of 
k • The procedure is repeated as 


n 
is increased until the desired period of time has been covered by the 


integration. 


The initial values of 
Q 
and 
U 
are chosen at 
tn 
rather than 


tn- 1/2 
in the above procedure, but little error is introduced since the 


initial values may be adjusted accordingly and the time step 
I::, t 
may be 


chosen as very small on the initial cycle. 


, 


The actual input model consis ts of a set of 
R, U, T, and 
V 


values for each zone of the configuration. Initial values of 
Q, E, P, 


Q, and :1 are then found by use of Eqs. (A. 12), (A .14), (A. IS)' and 


(A. 16) before starting the first time cycle. 


The set of Eqs. (A.9) through (A. 12) depend only on the localized 


properties of the fluid and they can be advanced explicitly in space and 


time, subject to the limitations on the time increments according to the 


Courant criterion (see Chapter 2). 


The set of Eqs. (A. 13) through (A. 28) must be solved simultaneously 


for all of the 
L 
zones in the configuratiofl because of the linkage between 


distant zones caused by the rad,iative flux. Since the advanced values 
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n+l 
V k-1/2 
of the specific volumes are known by advancing Eq. (A. 11) 


explicitly, a set of trial temperatures 
n+l 
T i _ 1/2' i = 1, ••• 
L 
are 


estimated, substituted into the Eqs. (A.13) - 
(A.28), and a Newton- 


Raphson iteration scheme is then used to adjust the estimated values 


of temperature until Eq. (A. 13) is satisfied to a predetermined accuracy. 


The iteration is carried out by numerical methods. The equation 


of state derivatives, at fixed 
V, are carried out by raising and lowering 


the temperature 2% from its trial value, e.g.: 


= 
E(V, l.02T) - E(V, O.98T) 


l.02T - O.98T 


The derivatives of the radiative quantities are computed by a 


ripple zone technique in which the temperature of a single zone is raised 


2%, the set of Eqs. (A .17) - 
(A.28) is evaluated and the desired derivatives 


formed, and the displaced temperature is then returned to its undisplaced 


value. This procedure is repeated for each zone in the configuration until 


all the desired derivatives are available. In practice, radiative derivatives 


for zones more than 2 zones distant are small, so that only a 5 zone set (,f 


derivatives is carried. Neglect of the more distant derivatives does not 


constitute a neglect of radiative transfer between distant zones; their 


neglect only influences rate of convergence of the iteration procedure. 


For a system of 
L 
zones, the iteration procedure then results in 
L 


linear algebraic equations in 
L 
unknown temperature increments, each 


equation consisting of a set of terms involving 5 of the unknowns. This 


array is then solved by direct elimination and back substitution. The 


temperature increments are then used to adjust the trial values of 


94 


(A. 29) 


temperature until all tempera~re increments for all zones are simul- 


taneously less than 10% of their respective zone temperatures. 


The FIREBALL code is written in a programming language called 


FORTRAN; the particular v-intage is known as FORTRAN II, Version 2. 


The various types of FORTRAN currently in use differ from each other 


only in minor details. This particular version was selected primarily 


because it has been in use for a number of years, and has achieved a 


measure of stability and reliability not to be found in the more recent 


efforts of the computer industry. 


The FORTRAN language little resembles the machine language coding 


of a decade ago, and its resemblance to ordinary algebra is so close that 


the average physicist or engineer can learn to read FORTRAN with a very 


minimum of effort. This allows one to communicate the solution of a 


complicated theoretical calculation to his fellow scientists in complete 


detail and complete scientific honesty. 


Section A.2 is devoted to a brief discussion of how to read FORTRAN 


and is designed for the scientist who is not familiar with this type of 


language. The remainder of Appendix A is devoted to the scientific aspects 


of the FIREBALL code and is intended to be independent of the language details. 


Readers not interested in great detail may thus skip over Section A. 2, while 


those interested in such detail will find this section helpful in reading the 


code itself which is listed in full. The four digit numbers which appear in 


parenthesis throughout this appendix, e.g. (0136), refer to serial numbers 


(line numbers) in the code listing. 
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A.2 
Reading FORTRAN 


Variables and Constants 


Algebraic variables are represented by symbols as in ordinary algebra 
-38 
~38 
and may take on values from about ±. 1.7 x 10 
to :!.. 1.7 x 10 
, and 


zero. Arithmetic is accurate to 8 significant figures. 


Variable names must consist of 1 to 6 characters, the first of which 


must be a letter of the alphabet other than I, J, K, L, M, or N. 


Examples: 
YIELD, X, A47, FI, DIVFAR 


Integer Variables 


Integers are used in subscripts, counting, indices, and sometimes 


in exponents. They are represented by symbols of 1 to 6 characters and 


must begin with 
I, J, K, L, M, or N. 


Examples: 
N, J63, MCOUNT 


Subscripts 


Variables may have up to three subscripts, but no superscripts. The 


su bs cripts may be positive integers, bu t not zero. 


Examples: 
X(K-tl), Z{J,N), BB(I,J,K) 


Arithmetic Operations 


The symbols for standard arithmetic operations are little different 


from those used in ordinary algebra: 
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x + Y 
means add 
X and Y 


X - Y 
means subtract 
Y from X 


X * Y 
means multiply 
X by Y 


X 1 Y 
means divide 
X by Y 


Normal algebraic sign conventions are used. A double asterisk is used for 


exponentiation, i. e.: 


R ** 2 
means 
R2 


R ** BETA means 
R~ 


Q ** 0.5 
means 
0 1/2 


Operations are carried out in the preferential order of first exponentiation, 


then multiplication and division , and finally addition and subtraction. 


Equal Siqn 


The equal sign has a meaning slightly different from algebraic usage, 


e. g.: 


XYIELD = XYIELD 
+ 
(A ** 2) - 
B 


means carry out the algebraic operations to the right of the equal sign and 


store the result of this computation as the new value of 
XYIELD. Note 


that the previous value of 
XYIELD 
is first used, then destroyed and 


replaced. 
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Parenthesis 


Parenthesis may be used to group information, indicate subscripts, 


or indicate the argument of a function. 


Examples: 


(AT 8) *C means add 
A 
and 
8, then multiply by C 


X(K+ 1) means 
K+ 1 
is a subscript of the variable X 


SINF(X) means 
X 
is the argument of the sine function 


Usually parenthesis will be used for grouping; subscripts can be 


recognized because they are integers, i.e., 
I, J, K, etc.; functions 


can usually be recognized by the letter 
F 
and the lack of an arithmetic 


operation symbol. 


Library Functions 


Certain library functions may be called in by name in order to save 


programming labor. 


Examples: 


Program Flow 


EXPF(X) 
= 


SQRTF(Y) = 


LOGF(Z) = 


x 
e 


log Z 
e 


Statements are processed by the computer in order of occurrence 


unless other directions are provided. Formula numbers, which occur in 


the first 5 spaces to the left, are not required unless control is to be 
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switched to that particular formula. The most common control commands (lre~ 


GO TO 436 


IF(X(K) - 3.0) 450, 500, 750 


sends control to statement 436 


tells the computer to test 


the sign of the expression 
in parenthesis (X-3) and to 
transfer control to statement 


450 if negative, 500 if zero, 


750 if positive. 


The IF statement provides the only means by which the computer 


makes judgments. 


Subroutines 


Complete subprograms, designed to-accomplish a particular set of 


computations many times, can be used. For example~ 


CALL STATE 


appearing in a program will send the computer to subroutine STATE, the 


set of computations indicated in that subroutine will be carried out, and 


control will return to -the main program at the statement following CALL 


STATE. 


Repetitive Operations 


To save labor in programming, the DO statement may be used: 


DO 100 
K= I, 9 5 
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instructs the computer to carry Out all instructions following this statement 


until statement 100 is reached. Do this first with K equal to 1 ; repC'at 


this process, increasing K by unity each time, making the final pass wi th 


K equal to 95. 


This statement can be used to compute repeti tive formulae I or it 


can be used to select values from an array of numbers. For example~ 


10 


50 


70 


TESTV = 
1. 6 E+O 2 


DO 5 a 
K= I, 100 


2 
means set TESTV = 1.6 x 10 


I F( x( K) 
- 
TESTV) 


XCRI T 
= 
X( K) 


50 I 
50 I 
10 


N = K 


GO TO 70 


CONTINUE 


XCRI T 
= 
0.0 


YSTART = XCRI T + DELTAY 


sets XCRI T equ al to X( K) 


retain the value of K as the symbol N 


exit from the DO statement 


a dummy statement 


set XCRI T to zero 


the next step in the computation 


selects the firs t X value that is greater than 160 out of a list of 100 


values of X, starting at X( 1) I •••••• x(lOO). If no such value is found I 


then XCRI T is set equal to zero. 


The above material has been brief and oversimplified I but should 


enable those unfamiliar with FORTRAN to decipher formulae from a program 


listing and to follow the scheme of computation. Additional help can 


usually be obtained from programmers if additional questions of interpre- 


tation should arise. Books on FORTRAN are available I but none are 


recommended. 
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A.3 
The computational scheme 


The general computational scheme is cyclic in nature, each completed 


cycle representing an advance in the time variable. 


The method of computation is illustrated in the schematic flow chart, 


Fig. A-2. Each block in this illustration represents a principal feature 


of the computation; numbers in the upper left corner of each block refer 


to a line number to the far-right of each page in the code listing. 


The computation begins by reading in an initial configuration from 


data cards. 
The configuration consists essentially of the initial density, 


temperature, velocity, and dimensions of each of the 100 zones in the 


model. The data cards may represent an initial configuration I or the data 


might represent the results of a previous computer run on a mOdel that is 


being done in short segments. 


If the entry is a restart rather than an initial start, a dummy subrc.utine 


rejust (0693)' is provided to make minor changes in program control 


parameters without the need for recompiling the major subroutines. 


After reading in the initial configuration, the data is tested for 


obvious errors (0367), the equation of state subroutine (038"7) I the radiati,,'e 


properties subroutine (0389) I and finally the flux subroutine (0410) are 


used to complete the details of the initial configuration. The computation 


cycle proper is then initiated (0402-0405). The initial time, which has been 


modified (0404) I is now restored to its proper value (0123). 


A sequence is then used to determine the number of zones, out of 


the possible 100 I that are to be used at the present stage of the computation. 


(0124-0138). The number of zones in use is continuously adjusted during 


the computation to avoid unneces sary computation. When the data is loaded, 
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zone 100 must have the ambient temperature and density cn burst ultitud8. 


The code then probes inward to a point where the temperature is tw.ic:,-~ 


the ambient value or the particle velocity reaches 1 m per second. The 


number of zones in use is then arrived at by adding a 6 zone safety factol 


to the result of the above selection. 


The next sequence (0139-0162) selects the times at which the datu 


print-out routine (0833) is called. selects the magnetic tapes to be used 


for the data. and sets up the necessary parameters to integrate the tOlal 


energy lost by radiation between two data printouts. and to find the radiant- 


power by differentiation of the energy vs time. The data print-out routine 


will be described in detail in Section A. 4. The data printout routines 


are one-way streets so far as the main stream of the computation is 


concerned. Data is siphoned ofL but no feed back into the computation 


proper occurs. 


The sequence (0164-0177) provides an emergency data saving 


mechanism for use in case the computation becomes numerically unstable. 


In order to provide an economical computation. the time step used must 


be just below the critical valuE.; prescribed by the various criteria which 


will be described in the hydrodynamic routine (0360) ~ These criteria are 


not fool"'-proof. and should the computation become numerically un stable 


(or should the compu ter operator err). the already completed compu toti Oil 


might be lost unless there were a mechanism for restarting. Once the 


calculation has become unstable. the data presently in the compu ter may 


be invalid. The present sequence writes the fireball configuration on a 


magnetic tape at the completion of each 50 code cycles. Each time another 


50 cycles are completed. the tape is rewound and the next configur3.tion 


102 


written. Finally, at the end of the run, the last configuration remains 


on the tape, is read into the computer (0305), and data cards are punched 


(0312). Should this data card punching process fail, the tape itself is 


saved and the cards can then be punched directly. At the time of terminati,m 


due to any type of error, data cards are then available not more than 50 


cycles back. On occasion, an instability may occur just after a configuration 


has been written in the tape, in which case the computation is lost. 
But 


this has a probability of occurrence of about 1 chance in 50. 


The next two sequences (0179-0257) select the rezoning subroutine 


(0179) and the zone splitting subroutine (0225). These sequences will be 


discussed below in Section A.5 dealing with these subroutines. In 


principle, one tries to remove fine zoning where it is not needed, and to 


create fine zoning when the physics of the problem so demands. Use of 


fine zoning throughout is not feasible because of the increased computational 


costs and also because the computer can only accommodate a total of 100 


zones at one time. If the finest zoning in the problem were used through Jut, 


the entire fireball could not be accommodated. So far as is known, other 


codes (Brode and Whittaker, private communications, 1965) currently 


in use are rezoned manually by visual inspection. The present sequences 


are an attempt to carry out this process automatically during the computation. 


The sequence (0259-0264) provides for the printing of detailed 


diagnostics at the initial cycle, the next cycle following, and at one 


selected by data card input. The diagnostic routine (1125) is called for 


this purpose. 
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A termination sequence is provided (0266-0342) which is used whenever 


the computation is terminated I except in those cases where the normal 


course of the computation is interupted by a machine difficulty or an 


operator error. The termination sequence provides for a final data printout 


(0274) I the printing of program diagnostics (0273) I the punching of current 


data cards (0275-0282) for possible restart purposes I completes I copies I 


combines I and unloads various tapes (0283-0326) I and punches data 


cards containing the configuratl"ons from 1 to 50 cycles prior to the 


termination (0303-0319) I before completing the run (0343). 


The series (0123) through (0345) has dealt primarily with control 


mechanisms I rather than the actual computation I which is resumed at (0346). 


In the course of the computation it is necessary to carry most of the program 


variables for all mass zones I but only for the current pOint in time. A few 


of the variables must be carried for all space zones I but for two consecutive 


pOints in time. For example I the present value of the specific volume for 


zone K is represented by the symbol V(K) I while the specific volume fromJne 


time- step in arrears I the retarded value I is represented by VR(K). As the; 


program cycles I the present value becomes the new retarded value I and a 


new "present" value is computed. 


The shifting process I in which the retarded values are set equal to 


the present values I takes place in the next sequence of commands (0346-0356). 


Also in this sequence I the present values of the variables are independeni.1y 


saved under the symbol W(K,I) for possible use in the event that the entire 


cycle requires restarting. The restarting procedure is discus sed at the end 


of this section. 
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Following retardation of the variables, a hydrodynamics routine 


(0360) 1s called, which will be des cribed in detail in Section A. 6. Thi s 


routine advances the velocities I radii, specific volumes I and artificial 


viscosities. 


One of the more important steps in the program is accomplished by 


a single dummy statement (0364). The first step in solution of the energy 


equation (Eq. (A .13», is the estimation of a new trial temperature for 


each zone. Skillful selection of the trial values will make the iteration 


process converge faster and thus speed up the computation. However, in 


an iteration process the equation being solved is never satisfied identical:y I 


and any prejudice used in estimating the new temperatures may tend to 


impose non-random errors in the final results. For this reason, the old 


values of temperature are used as the first estimates for the new values. 


The iteration scheme itself then, in a sense, becomes a basic part of the 


system of equations. As will be discussed later, the iteration scheme used 


employs the same set of equations as are used to represent the energy 


equation itself. This procedure should tend to minimize any bias. 


The iteration cycle proper (0365) begins by testing the set of 


temperatures for negative or zero values. Such values are likely to occnr 


in the event of a numerical instability I in which case the computation is 


immediately terminated through" the sequence previously described. 


In the next sequence (0371-0387) the equation of state subroutine (1217) 


is called on three successive Qccasions. This subroutine is discussed in 


Section A. 7. At this pOint in the computation, the advanced values of 


the specific volume are known for each zone. A trial set of temperatures 


has been selected. The derivatives of the internal energy and pressure, 
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with respect to temperature, at constant specific volume, are formed by 


first increasing the temperatures 2% above the trial values, passing through 


the equation of state to obtain energies and pressures, then repeating the process 


after lowering the temperatures 2% below the trial values. The required 


derivatives are then formed numerically from the above data (0384-0385). 


The temperatures are then returned to the original trial values and a 


final pass through the equation of state is made to obtain trial values of 


the internal energy and pressure (0387). 


Two subroutines are then ':called to provide the neces sary radiative 


flux divergences for use in the energy equation. Subroutine SWABZ (0490) 


uses the values of V, T, and the dimensions of the zones to compute 


values of the SWAB and Z functions and is discussed in Section A.8. These 


values are used in subroutine FLUXS (1183) to compute fluxes and flux 


divergences as described in Section A. 9. 


A temperature test occurs next in the computation, but is bypassed 


on the initial trial of the energy equation. 


The energy equation itself is then evaluated (0415-0424). A residue, 


comprising the imbalance of the energy equation due to the use of trial 


temperatures, is computed for each zone (0422). 


The temperature iteration of the energy equation then commences ' 


(0425). Two types of iteration schemes can be selected. If the configuration 


is entirely optically thin (i. e. the total optical depth from the outer edge 


to the center is less than a given value) a non-radiative iteration (0429) 


can be used. Radiative transport is still included in the energy balance, 


but explicit radiative derivatives do not appear in the iteration scheme. If 


this alternative is selected, the temperature increments for each zone are 
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computed immediately (0431). 


When radiation is included in the iteration (0439), the radiative and 


hydrodynamic derivatives are ccrried only in the inner zones, while only 


hydrodynamic deriVatives are carried in the outer zones. 


In the radiative iteration procedure, subroutine COEFF (1395) is 


called to form the required coefficients for the solution (0447) for the 


temperature increments. This subroutine is described in Section A. 10. 


Once the full s·et of proposed temperature increments is known, they 


are tested (0457) to see if they are within arbitrary bounds which have 


been developed by e1tperience. Should these bounds be exceeded, there 


is a high probability that the time increment used by the computer was too 


large. Program diagnostics are then printed, the variables are restored 


to their vall1es at the beginning of the cycle, and the cycle is repeated 


using a smaller time increment. This recycling is allowed 3 times, after 


\ 
which the computation proceeds even though the test bounds were violated. 


If the test criteria were correct, a numerical instability results, and the 


computation is terminated from ·:me of the several pOints in the program 


where the instability can be positively detected. An instability does 


not always occur, however I since the tests are not infallible. 


Usually a recycling occurs when a time increment just slightly too 


large has been used and a single recycle cures the difficulty. Use of this 


stability check "after the fact" enables the various time step selection 


criteria to be pu shed close to their limits. Use of large safety factors 


in these criteria, as is the common practice, would be prohibitive as 


the total computation time woul? be seriously increased. If the proposed 


temperature increments satisfy the stability tests, they are accepted and 
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the trial temperatures are modified accordingly (0485). 


These new temperatures become a new set of trial temperatures 


and the iteration cycle is begun again (0367). The iteration convergence 


test (0409) is made well 1nto the second and subsequent passes through 


the iteration sequence. Should the test be satisfied, then the best values 


of internal energy, pressure, ~d fluxes, etc., are available as the cycle 


is completed. If the test is failed, the iteration cycle then progresses with 


all the necessary data. The code as written performs some unnecessary 


computations in that the state derivatives are computed but not used if 


the convergence test is satisfied. This technique was a compromise 


between several alternatives and was adopted because it required less 


computer storage at a period when the program was storage limited. 


The actual convergence test used (0408) consisted of the requirement 


that the temperature increments for all zones be simultaneously less than 


10% of their respective zone temperatures. This method proved to be 


more effkient than several other methods tried which were based directly 


on the degree of imbalance of the energy equation. 


The number of passes through the iteration sequence is limited to 


3 on anyone code cycle (0416). It is a characteristic of the Newton method 


that convergence takes place very rapidly or not at all, and that on odd 


number of attempts is usually better than an even number. Should an 


abnormally large temperature increment occur after 3 iterations, a 


convergence check at the beginning of subroutine HYDRO detects this 


and forces a smaller time step on the following cycle. This test differs 


from the main program stability check in that it is applied only after 


3 full iterations and uses stricter criteria. 
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A.4 
Data printout routines 


The data printout sequence is contained in subroutine CGSPO (0833), 


This routine samples the data from the computation proper, generates 


additional parameters of interest from this information, and produces tW(· 


magnetic tapes: the main listing and the user tape. 


The general problem of determining observables from a list of 


temperature and density values is a difficult one. Many of the "observable" 


parameters generated by the printout routine are extremely crude and must 


be used with extreme caution. These "observables" have been printed 


below the main tables of data and include the shock radius, fireball 


radius I effective temperature, color temperature, and the spectral 


distributions. For example, the code defines the fireball radius (FBR) 


as the radius at which the optical depth is 0.44 as measured from outside 


the fireball. This simple definition is easy to code and is useful radius 


to have printed out. att for comparison with an experimental radius 


measured photographically, a complete brightness profile is required for 


comparison with the corresponding densitometer traces from the experiment. 


The routine begins by testing for temperature inversions (0864) and 


causes program diagnostics to be printed 1£ such inversions are found. 


Placement of this test in the printout routine causes the diagnostics to 


coincide in time with the data printout so that any unusual features seen in the 


listings may be studied in detail. 


In the sequence (0868-0874) a number of variables to be printed out 


are set up. This sequence works in conjunction with the sequence (0139- 


0162) in the main proqram. The variable CN(LZ) is the average time step 


between two succeillive printouts. The variable SAVEl is the value of 


the last time step before the printout. The total power being radiated 
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in watts (POWER), the color temperature (TGOL), and the effective 


temperature (TEFF) all depend upon the value of the total flux (TFLUX) 


which is the power being radiated in ergs. TFLUX is computed in the 


main program (0155). The instantaneous total flux at the end of each 


cycle is called FLOX and is summed continuously (0139) throughout the 


computation and called FLEX. The average power TFLUX over the time 


interval between printouts is obtained by numerical differentiation of 


FLEX (0155). 
This method gives more representative results for the 


power output since the instantaneous power tends to fluctuate due to the 


finite zoning of the model. The instantaneous total power corresponding 


to the time of a given printout can easily be obtained since the radius 


and outward flux at the last zone are available from the printout. 


The code definition of the effective temperature (TEFF) is defined 


(0874) from 


TFL UX = 4TT (F BR) 2 a (T EFr) 4 


where 
a = 
-5 
5.67 x 10 
= Stephan's Constant. 


It is the temperature of a" black body radiation having the same size 


as the fireball and which emits the same total power". This temperature 


is a minimum estimate of the fireball "surface" temperature sin ce the 


fireball does not have an emissivity of unity. 


The code definition of (0873) color temperature (TGOL) uses an 


estimate of the total emissivity and yields a higher temperature which 


can be used to describe the sp~ctral character of the radiation escaping 


+ 
from the fireball. The emissivity estimate is based on the Z 
function 
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for the zone having the largest ')., 
and the radiation temperature for 
c 
that same zone (0689). 


After the fireball becomes transparent in the continuum, the effective 


and color temperatures are no longer defined since their definitions 


involve the fireball radius, and by the code's definition, FBR=O under 


these circumstances. 


At this pOint (0875), the printout routine calls subroutine PHOTOG 


(1061). This routine computes an estimate of the photographic brightness 


as a function of radius. Using an approximate fit (1073) to the absorption 


coefficient for the "photographic II region of the spectrum, and a Planck 


function (1083), the emitted intensity is calculated for chord rays which 


are tangent to the mid-point of each mass zone. Since the spectral band 


width of the, photographic region is unspecified, only relative brightness 


values are obtained. The brightness scale, however, remains fixed 


throughout the entire computation. It is clear that only an estimate can 


be obtained from this comput8'tj,on since the absorption coefficient is not 


independent of wavelength across the photographic region of the spectrum, 


and also because the appropriate absorption coefficient data are not 


available at temperatures above 20 eV. The data is useful, however, in 


making comparisons Wi th exper!mental results and in finding the gross 


brightness variations across the fireball. 


, 


The CGSPO routine next makes an estimate of the spectral distribution 


of the total radiation emitted by the fireball. The code decides (0877) 


whether the fireball is optically thick or is transparent. If the fireball 


is transparent (0878), no realistic estimate of the spectral distribution 


can be made since the spectrum consists of emission lines and a weak 
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continuum. An extremely crude estimate of the "visible power" 


(4000-7000 A) can be made by noting that fireballs tend to form an 


isothermal region at these late times. Most of the emi tted radiation 


must come from this isothermal region since it is the hottest part of 


the fireball. These isothermal region temperatures range from 9500 0 K 


down to 5000 0 K as the transparent fireball cools. If it is assumed that 


the envelope of the emitted spectrum crudely approximates a Planckian 


distribution, then 37 -t 3% of tha radiation will fall in the "visible" region 


of the spectrum over this entire range of temperatures. Hence, for the 


transparent case, the visible power (P47) is taken as 37% of the total 


power (0878). 


If the fireball is optically thick (0888), then the color temperature is 


used as a basis for division of the total power into broad spectral bands 
.. 
(0889-0939). The formula for the Z function, which very nearly approximates 


the fractional Planck function in the visible and IR spectrum, was used in 


<Xl 


place of an accurate fit to J ~ d). 
as an analytic fit to this integral was 


not know to be readily av~ilable. The accuracy so achieved· i s probably 


better than the physics of the spectral estimation proces s. The overall 


results are of a quality comparable to those achieved by assuming the 


sun to be a 6000 0 K black body. 


The next sequence (0944-0963) performs an energy check by summing 


the internal and kinetic energy of the current configuration. The sequence 


computes the internal energy and the kinetic energy zone by zone, and 


also the internal energy of the volume of space currently occupied by the 


configuration if there had been no detonation. This "ambient" energy of 
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the undisturbed air constitutes a significant part of the total energy of 


the configuration, particularly at late times when the radius of the 


configuration is large. It must be taken into account in studying the 


energetics of the model. 


The code has no built-in method of forcing energy conservation. 


The energy check sequence is simply an after-the-fact sampling to see 


if energy has been conserved. 
If the difference equations do in fact 


represent the basic conservation equations, then mass, energy, and 


momentum will automatically be conserved. 


The energy accounting at anyone moment can be verified by adding 


the present model energy due to the detonation and the thermal losses up to 


the given moment. The quantities TOT ES, EAMB and TYIELD are printed 


out and are the appropriate numbers to use for this purpose. TOT ES in 


the printou t i~ the total energy, (DPP(LZ) in the code (0959)) and includes 


EAMB, the energy content of the cold air before the detonation. TYIELD 


in the printout, ~LEX in the code (0139)) is the total energy radiated 


outward across the outer boundary of the model from the time of 


detonation up to the present time. In general, the code conserves energy 


to within a fraction of 1%. 


Sandwiched in the energy check routine are some operations which 


concern the printing out of shock parameters. (0955-0958). It should be 


remarked that the location of shock fronts in a numerical configuration is 


a difficult problem because of the varying number of shocks and the wide 


range in their characteristics. Many methods were tried, but with limited 


succes s. Shocks are best located by a careful study of density vs. radiu s 
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and velocity vs. radius graphs, 


The location of a strong shock is determined by the program by 


finding the largest value of the artificial viscous pressure, 
with the 


added requirement that it be larger than the ambient gas pressure at 


burst altitude. This shock radius is printed, but cannot be relied upon 


without further inspection to see if a shock really exists at that radius. 


This simple technique will print out the location of one shock, but it 


may not select the same shock. on consecutive printouts. 


The sequence (0978-0991) 'imply changes units and sets up certain 


quantities in proper form for priMing them out. For example, (0987) 


computes the effective value of "gamma, " while (0984) yields the 


temperature in electron volts. 


The sequence (1001-1016) prints out the main listing, while the 


sequence (1022-1029) prints out the user magnetic tape for use in 


continuation computations such as fireball environmental studies. 


A.5 
Variation of the number of zones in use 


It has been said that there are three major problems to be solved in 


writing a fireball code: the basic physics, the interpretation of the 


results, and the use of proper zoning. Proper zoning is by far the most 


difficult of these problems. One simply cannot use a fine mesh throughout 


the configuration as the computing time increases approximately as the 


square of the number of zones, i.e. halving a zone halves the time step 


and doubles the number of computer operations I thereby requiring four 


times as much computer time. 
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In the other extreme, if the model is represented by zones of equal 


size, some of the zones may be too large to properly represent the 


gradients in phys ical properties". 


Consider, for example, a model where the fireball has a radius of 


103 cm and may require zones of width 10- 1 cm at the fireball boundary. 


Since the computer can. accommodate only about 100 zones at a time, the 
-1 
zone obviously cannot be of equal size. Furthermore, use of a 10 
cm 


zone in the hot isothermal region would slow the computation by at least 


a factor of a hundred, and perhaps a factor of ten thousand or more. 


Ideally, new zones should be created and old zones combined in an 


optimum manner based on the local physical conditions at each pOint in 


the" configuration. The subroutines SPLIT (0785 
) and REZONE (070 1) 


are provided for these purposes. 


In the split subroutine, the zone indices for the zones exterior to the 


zone to be split are first shifted outward (0789-0808) to make room for 


the new zone. 


The zone is then split in half mass-wise (0810-0811). The particle 


velocity at the new interface is taken as the square root of the average of 


the squares of the particles velocity at the boundaries of the original zone. 


The temperature of the two new zones are displaced 10% above and below 


the temperature of the original zone, (0824-0825) and a similar treatment 


is accorded the internal energies (0826-0827). 


This routine must be considered as only a crude beginning and much 


work is still being done in developing new techniques. The present 


routine can only be used before the gradients become large, and even then 


its use should be discouraged. 
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The calling sequlnce f« the SPLIT subroutine occurs in the main 


program (0225-0257). Statements (0228 and 0248) prevent a split more 


often than each third eycle. Statement (0231) requires a split at any time 


that a single zone represents more than 8% of the radius of the entire 


conU9\lration. 


The real purpose of zone splitting is to provide very fine zones just 


ahead of an advancing shock front so that the optical properties will not be 


too severely distorted. The sequence (0232-0244) attempts to achieve 


this goal. Optical depth 0.7 from the outside is located (0232) and this zone 


is tested to see if it is being compres sed (0 23 5). Then if the particle 


velocity is greater than 105 cm/sec as is characteristic of optically thick 


shocks, the five zones immediately exterior are scanned and one per cycl8 


can be split, until a certain minimum zone size is reached. Once the 


shock. has started I the splitting takes place 5 zones ahead of the shock 


until til. shock becomes transparent. The minimum zone size has not yet 


been formulated in general terms I and must be re-programmed in for eaC1 


separate run. In the program listing (0240) a zone of 200 cm or more can 


be split, resulting in a minimum zone size of 100 cm. This is about 


appropriate for a megaton burst at sea level. 


As a final insurance against undesirable zone splitting at high altitudes I 


statement (0244) prevents splitting at temperatures above 6 x 10 4 oK. 


Subroutine REZONE (0701) I which combines two existing zones is on 


a much more sound basis. Conservation of mas s in the zone combining 


process is trivial (0735) I but the conservation of energy and momentum are 


slightly more complicated. The internal energy of the new zone is taken 
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as a mass weighted average of ti1e internal energies of the two zones being 


combined (0733). The conservatism of kinetic energy and momentum for a 


four zone system (0709-0712), which is collapsed into a 3 zone system, 


results in two equations for the two unknown particle velocities at the 


inner and outer edges of the combined zone. 


The remainder of the subroutine REZONE consists of a shifting of 


indices for zones exterior to the fusion, and the addition of a new zone 


100 at the outer edge of the configuration. 


While this routine works well and quite accurately. some skill is 


required (but not always attained) in deciding when rezoning should take 


place. 


The calling sequence for rezoning is in the main program (0179-0224). 


On the basis of optical properties (0182), the outer limit of the region to 


be seaMed (the index NZ TS) and an allowable mass ratio for neighboring 


zones (PMT) are selected (0181-0202). If at least two cycles have passed 


since the la.st rezone (0204), zones 9 through NZTS are tested for size 


compared to the size of the entire configuration (0205). temperature 


gradient (0206), density gradient (0207), and mass gradient (0208). If 


these gradients are less than the allowable artitrary limits. then rezoning 


is allowed so long as the zone is not undergoing significant compression 


(0209) • 


A.6 
Hydrodynamics routine 


The hydrodynamics routine (1261) is patterned after the artificial 


viscosity method of treating shocks (see Chapter 2~). In addition to 


advanaing four of the basic equations, this routine also controls the size 


of the time increments and performs miscellaneous other functions. 
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At the beginning of the routine I the previous temperature increment 


is tested (1273) against the temperature, and if the increment is too large, 


the proposed time step is cut a factor of five (1277). This test serves 


several purposes depending upon the values of the previous temperature 


increment that may be stored in the computer at the moment. If a mainc'ycle 


has been completed' satisfactorily, 1. e. all temperature increments are les s 


than 10% of their respective temperatures, then no decrease in estimated 


time step takes place. If, however, a main cycle was completed after 


three passes of the energy equation, the 10% requirement may not have been 


satisfied, 1. e., the iteration may not have converged. Should this be the 


case, then the time step is shortened in the hope that a numerical 


instability can be avoided. 


If the main program stability check (0457) is violated, the program 


returns to the beginning of that master cycle after restoring all data 


except the temperature increments. Then, since the HYDRO tests are 


slightly more demanding than are the main program tests, the time step 


will also be decreased when the main program senses a difficulty. 


The equation of motion sequence is carried out in two separate 


phases which are separated by criteria to choose a 'new time increment. 


The particle velocities are first advanced using the previous value of 


the time step (1289) so that proper time centering of the difference 


equation (Eq. A .9» is maintained. 


A new time increment is then selected (1298-1335). First a time 


step 30% larger is suggested (1299). This value is then reduced should 


the Courant criterion, Eq. (2.5-4), demand that a smaller value be used 
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(1302-1310). The time may be further reduced by a radiative criterion 


(1314-13 26). Gross checks are then imposed, as a safety factor, which 


demand that the new time step can never exceed the previous one by more 


than a factor of two, or be less than a given minimum value (DTMIN). 


The value of DTMIN is continually increased during the computation and 


is kept a factor of 50 smaller than the largest time step that has been u~ ed. 


Should the speed of computation decrease more than this factor of 50, too 


large a time step (the value of DTMIN) is forced into use and the computation 


may become unstable and turn itself off. This feature prevents the waste 


of computer time should the time step conditions become abnormally 


critical at some pOint in the configuration. 


Superimposed on these criteria is the mechanism for causing data 


to be printed out exactly at fixed predetermined times (1333-1335). This 


criterion may further decrease the time step. 


After a new time step has been decided upon, new values of the 


particle velocity are determined by linear interpolation (or extrapolation) 


so that proper time centering of the equation of motion is maintained (1337). 


The radii (1344) and specific volumes (1362) are then advanced in a 


straightforward manner. The new radii are tested before adoption (1348) 


to prevent sudden zone collapse should the estimated hydrodynamic time 


step be too large. A local recycle with decreased time step is then 


instituted (1358). 


It is perhaps appropriate to again mention that the use of large 


safety factors in the time step criteria result in the use of large amounts 


of computer time. Considerable economy is achieved by lowering the 


safety factors, testing the proposed results before closing a cycle I 
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and re-cycling when necessary. Using this technique and with a little 


experience I the code can be kept running at near optimum speed. 


The sequence (1372-1385) advances the artificial viscous pressurE. 


A variety of formulae are available in the literature for this purpose. The 


formulae used here is similar to that given by Richtmyer (see Chapter 8 


references). This form appears to give better results when shocks reflect 


at the center of the sphere than do the "linear" and "quadratic" forms. In 


the final analysis I the use of an artificial viscosity is an art and an 


adjustable constant (1374) is available to achieve optimum results for 


any specific situation. Use of too small a constant causes numerical ripples 


behind the shock I while use of too large a value spreads the shock-over too 


large a distance and lowers the shock velocity. According to Richtmyer I 


choice of the arbitrary constant in a manner so as to spread the shock 


discontinuity over 3-4 zones results in shock velocities and pressures 


that agree well with laboratory experiments. 


The use of an artificial viscosity causes a ficticious precursor ahead 


of the shock front. If this precursor is optically thick I the shock radiation 


rate will be affected. This effect reduces the rate of radiation los s when 


the shock temperature is large I but numerical tes tirig has shown that it 


has little effect on the shock energetics. 


A.7 
Equation of state routine 


The equation of state subroutine (1217) is entered with known values 


of the temperatures and specific volumes for each zone I and values of 


pressure and interval energy are computed. 
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The source data for the analytic fit to equation of state are 


primarily due to Gilmore (1955) and to Hilsenrath and Beckett (1955). This 


polynomial type of fit, though in principle not as accurate as an iterative 


routine, is preferred for computational purposes because the derivatives 


are well behaved. This fit is also self consistent in that the hugoniots 


are closely satisfied, while some of the piece-wise fits, that are accurate 


over limited ranges, fail in this regard. No significant inaccuracies 


are known to have resulted from use of this simple expression. 


There is some question as to whether radiation pressure and the 


radiation energy density should be included in the equations of state 


(see Chapter 2). The subroutine STATE allows an explicit choice to be 


made in this regard (1222-3). These effects can only be important at high 


temperatures. The argument against including these effects is that at the 


high temperatures the matter and the radiation cannot be in equilibrium 


according to the local high temperature since the g as is too tran sparent. 


The radiation field is characteristic of the lower temperatures where the 


opacity is higher. To base the radiation pressure on the higher values 


of temperatures in the interior would thus be an overestimate. 


A .8 
Radiative properties routine 


The radiative properties of the fluid, which are characterized 


bytheS, W, A, andZfunctions , Eqs. (A.17) through (A.2S)' are supplied 


by subroutine SWABZ (0492-0692). The analytic representations of the Z- 


and A functions I which differ only in that the electron scattering is omitted 


in the A function, are formulated in terms of AH I a radiative mean free path 


which depends upon the zone temperature, zone specific volume, and the 


time. 
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The Z 
function (ZI in the code) is formulated directly as a 
* 
transmission (0500-0507). 
The A function cannot be formulated directly 


in this fashion since In A is needed to define optical distance, and 


truncation errors would occur as A ..... 1.0. Thus 
AH' 
(HMFP in the code, 


is found (0508-0511), together with a separate high temperature emission 


term (0512) due to Bremsstrahlung (HMFF in the code). The smaller of 


these two mean free paths is used for temperatures above 10 eV (0514). 


The emission optical depth (0516) then follows immediately from the 


zone thickness. The A function (0517) at this point in the code is identical 


wi th the definitions used in Eqs. (A. 21) and (2 .5-14). A later sequence 


(0670-0675) redefines A as (I-A) (an emissivity) as a computational 


convenience for zones of small emissivity. 


The variables Be and BS 91518 and 0526) are used to represent the 


source function at the zone centers and zone boundaries respectively. It 


will be noted that the boundary value is determined by linear interpolatior 


of the source function in terms of geometrical rather than optical dis tance. 


The two methods yield substantially the same results when the two neighboring 


zones have comparable optical thicknes ses, but the geometrical method 


gives better values when an optically thick zone occurs next to an optically 


thin zone. The terms involvinq th e boundary values cancel if both zones 


are optically thick, while the source function gradient is small for two 


transparent zones; SO that the method of interpolation is not important in 


these cases. 


* 


The sequence (0536-0563) computes the optical depth inward from 


The double subscripting, 1. e. A(K, N) has been carried over from an 
earlier version of the code in which N spectral bands were used. 
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the outer boundary of the configuration and sets various control variables 


used by the program at other locations. The index LZR specifies the las t 


zone to be included in the radiative part of the iteration scheme. The 


fireball radius, FBR, is set as the outer radius of the zone in which the 


total optical depth as measured from the outside ha s reached 0.44. ChOice 


of this value has no particular significance. For high air densities the 


optical depth increases abruptly so that any reasonable test value will. 


result in choosing of the same zone to specify the radius. For rarefied 


atmospheres, considerable limb darkening takes place and chord integrations 


would be necessary to find a precise radius. The test value of 0.44 has 


been found to yield radii satisfactory for the intended applications elsewhere 


in the code. The index MCP is used to specify the zone in which optical 


depth unity occurs as measured along the representative ray. (This 


corresponds to a radial optical depth of 2/3). 


The sequence (0564-0607) computes the wavelength at which the 


principal spectral absorption edge occurs for each zone. This wavelengtt: 


depends only upon the physical characteristics of each zone, i. e., its 


temperature, density, and zone thickness. Zones are first sorted 


according to temperature to determine whether the spectral transition 


is due to an atomic species, the nitrogen molecule, or the oxygen molecule. 


The transition edges due to atomic species are in the far ultraviolet 


and a fi t will be required only at very high altitudes when the fireball is 


transparent throughout while the temperatures are high enough so that the 


Planck function is significant tn the ultraviolet. This fit used by the coqe 


is therefore very limited in its application. 
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The nitrogen and oxygen molecular formulations involve a continuum 


with a relatively sharp long wavelength limit, and an extensive system of 


molecular bands extending from this limit toward longer wavelengths. The 


continua 
are assumed to exist at all times that the molecules have 


sufficient population, but the band systems, which depend upon the 


population of the higher vibrational levels, are assumed not to exist until. 


the pas s age of a vibrational relaxation period. 


While in principle the relaxation time should be measured from the 


time that each zone is first heated (actually the temperature-time history 


should be taken into account), a satisfactory first approximation is given 


by use of the time since the detonation. The relaxation time for N2 at 
-7 
-7 
sea level is taken as 5 x 10 
sec and for 02 as 3 x 10 
secs. These 


times are scaled inversely with the density to obtain the relaxation 


times at higher altitudes. These times are based on data by Blackman (1956). 


The absorption edge and the continuum absorption due to a particular 


species must fade away as the population is diminished by dissociation. 


At a given density, the population decreases very rapidly with temperature 


as soon as kT becomes comparable with the dissociation energy. Due to 


the finite zoning structure, the population of a particular specie will be 


appreciable in one zone, but negligible in a neighboring hotter zone due 


to the rapid temperature dependence of the populations. Thus, for the 


purpose of computing absorption edges, the species can be assumed to 


exist below a certain temperature and not to exist above that temperature. 


The dissOciation is thus assumed to take place at a temperature rather than 


over a narrow temperature range. The dissociation temperatures for N 2 
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and 02 ' as a function of density, are represented in the code by two 


formulae (1708-1709) and are used in the SWAF5l subroutine (0572-0580). 


+ 
The sequence (0612-0669} computes the Z values (ZO in the code). 


The Z+ function depends on 
A c 
which has been determined above, and 


the temperature 
TR 
used to describe the spectral distribution of the 


radiation incident upon the zone. The sequence (0622-0642) selects this 


temperature for each zone. In the case of an optically thick fireball, this 


temperature is that of the first zone at optical depth greater than 0.7 from 


the zone in question (0628). For a transparent fireball, where the radiation 


comes from a shell rather than from a "surface", the temperature is taken 


as an optical depth weighted average of the zones interior to the zone in 


question (0624). This "shell source" sequence is important only at high 


altitudes where the fireball is transparent, and when the temperatures are 


still high enough to have radiative flow in the far UV region of the spectrum 


where the absorption edges occur. 


+ 
The Z 
values are corrected for the effects of intervening zones 


between the source of the radiation and the zone in question (0655). 


The sequence (0676-0688) becomes effective at very high altitudes 


when the configuration is quite thin and thermal radiation plays only a 


minor part in determining the temperature distribution. Under these 


circumstances, the series expansion of the source function should not 'be 


truncated and derivatives of higher order than the first should be included. 


But since the radiation is of little importance in this case, the series can 


be further truncated SO that only the zero order term is used and the code 


automatically switches over to radiative transfer for isothermal optically 


thin slabs. 
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The variable EMS (0689) is used as an estimate of the total emissivity. 


+ 
It corresponds to the Z function for that zone having the maximum cut-off 


wavelength, but without correction for intermediate zones. 


A.9 
Radiative flux routine 


Since all the needed data has been generated elsewhere by the code, 


the radiative flux routine is particularly simple (1183). 


To perform the flux integration a boundary condition is needed. For 


this purpose, the inward flux incident on the outer boundary of the 


configuration is chosen as zero (1190). 


The inward flux at the boundary, one zone inside, is given by the sum 


of the transmitted and emitted components (1193) as in Eq. (A. 27). The 


integration continues inward until the spherical central zone is reached. 


The outward flux for this zone is the sum of the transmitted inward flux 


and the local emission (1195). The integration then proceeds outward (1197) 


using Eq. (A. 26). One can view this process as an integration starting3t 


one side of the fireball, passing through the center, and then on out the 


other side. All of the inward and ou tward fluxes needed to form the radial 


component of the flux divergences (1211) are generated during this integration. 


It should be noted that this subroutine is written again for N spectral 


bands, and that the total directional flux at a given boundary is obtained by 


a summation over the spectrum, even though only 1 spectral band is used 


by the present code. In addition, the A function used in (1193-1198) does 


not correspond to the transmission function used elsewhere in the text, 


but to (1-A) as explained in the SWABZ subroutine. 
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The energy Equation (Eq. ( A.5 » is solved by an iterative method as 


described by Hillendahl (1964). Subroutine COEFF (1395) performs the 


necessary algebra to form the required derivatives. 


The basic approach is to compute the derivatives numerically by a 


ripple zone technique, rather than to derive and to code explicit formulae 


for the derivatives. In using this technique, the temperature of a single 


zone is increased 2%, then formulae identical to those used in subroutine 


SWABZ are used to calculate those radiative properties which change as a 


result of this temperature increment (1429-1576). A flux integration is then 


carried out (1577-1593), the flux derivatives are formed (1614), and the 


perturbed variables are returned to their normal values (1627) which were 


saved (1411). This ripple zone process is repeated until all the required 


derivatives are available. Derivatives with respect to zone temperature 


more than two zones distant are truncated. 


The energy equation derivatives, including both the radiative and 


hydrodynamic parts, are formed numerically (1647). The data is then 


available to form a set of linear algebraic equations for the temperature 


increments of the zones. 


This matrix is then solved by direct elimination and back substitution 


using recursion formulae (1673) and (0447). 


This numerical method of solution, involving the ripple zone method 


of obtaining flux derivatives and the step by step formation of the energy 


equation derivatives has been found to be both convenient and economical. 


If details of the radiative properties or flux formulae are changed, 


duplicates of the new formulae are simply inserted into the coefficient 


subroutine without the necessity of deriving explicit formulae 
for the 
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derivatives I which in some cas8s I ,have as many as 40 terms. Such lengthy 


formulae are difficult to derive and code without errors I and the sorting of 


terms to avoid arithmetic truncation difficulties is a major task. Since the 


actual number of operations to be carried out by the computer is approximately 


equivalent in the two methods I the numerical technique is to be preferred. 


A.I0 Code listing 


The pages following presfmt a complete lis ting of an actual working 


radiation-hydrodynamic code which has been described in sections A. I-A. 9. 


No attempt has been made to edit the listing for publication purposes since 


this practice quite often results in the publication of codes that do not work. 


128 


I 


• 
FORTRAN 
0001 
CFIRE 1 
0002 
C 
HOOIFIEO OII1ENSION STATEMENT 
NSZ-l 
ONLY 
0003 
DIMENSION 
T(lOl). 
TS(lOO). 
TP(lOO). 
TI1(lOO). 
TR(lOO). 
0004 


2 
PAClOO). 
P(lOOl. 
PSCI001. 
PP(lOOl. 
PMClOO). 
PRCIOO). 
0005 


3 
OECIOO). 
EClOOl. 
GMCIOO). 
EPUOO). 
EI1ClOO). 
ERCI00). 
0006 
4 
Q(1OO).VCIOOl. 
VSCI00). 
RHOZCI00).UClOO). 
VRClOO). 
0007 
50RClOOl.OCl).RClOl). 
RZCI00). 
R2(100). 
UA(100). 
URCtOO). 
0008 


6 
OTM(100). DTRCIOO). OTHCI00). DPT(100). OET(100). RR(100). 
0009 


7 
ZI1AS(100).DIVFA(100l.DIVFRCIOO).RESDUE(100).FPelOO). FZetOO). 
0010 


8 
FI1CIOO). 
FMM(IOO). HZ(100). 
HP(100). 
HPPP(100).HPPCI00). 
OOtl 


9 
OI1"CIOO). 011(100). 
DZ(100). 
DP(100l. 
DPP(100). BN(100). 
0012 


1 CN(100).DN(1001.ENCI001.ETA(1001.XSCI001.GNUCI00).PART(100) 
0013 
DIMENSION THETACI00) 
0014 
C 
00t5 
C 
STORAGE FOR RADIATION VARIABl.ES 
0016 
C 
0017 
DII1ENSION 
88Cl00010). 
DTAUClOOolO).S(100.2). 
FOCtOO.lO). 
0018 
1FOHIOO.IO).FIC100.10). 
FH(100.10). 
ZO(100.10),ZI(100.10) 
0019 
2.AClOOolO). 
W(1000l0). 
9C(101010). 
FIS(100). 
FOS(100), 
0020 
3CWl(100. 2).YI1(100l,TAUCtOO. 2) 
0021 
C0l1110N T.TS.TP.TI1.TR.PA.P.PS.PP.PI1.PR.DE.E.GI1.EP.EI1.ER.Q.V.VS.RHO 
0022 
2Z.U.VR.DR.D.R.RZ.R2.UA.UR.DTI1.DTR.DTH.DPT.DET.RR,ZI1AS,DIVFA,DIVFR 
0023 
3.RESDUE,FP.FZ.FI1.FI1I1.HZ,HP,HPPP,HPP,DMI1.DI1,DZ.DP.DPP.BN.CN.DN.EN 
0024 


4.ETA.XS.GNU.PART.THETA.KGU.KGP 
0025 
5.BB.DTAU.S.FO.FOT.FI.FIT.ZO.ZI.A.W.BC.FIS.FOS,CWL.YI1,TAU,NSZ, 
0026 


6 NZSS.·. LZ.lZI11.lZI12.lZPl.LZP2,BR.RJ.RX.NCW.NI1C.NTI,FLOX,FlEX. 
0027 
7TII1E.DT.CS.CR.RS.I1Cl.I1CP.I1CW.RI1.VD.ZA,NB.NQS.FBR,LZR,DTI1IN,NR 
0028 
8.KZI.KZ2,KZ3.KZ4.KZ5.KZ6.KZ7.KZB.KZ.TD02.TDN2,NDP.SCAlE.EMS. 
0029 
9BlANK.AST.TEE.PLUS.PERIOD.DASH.EQUAL.PINUS.FFF,UUU.PPP.NTAPE.TII1EW 
0030 


6.TFLUX.FLIX.TIMES.P03.P34.P4S.PS7.P71.P47,Q47.Q71.WKT.YIElD.XYIElD 
0031 
KZ6=O 
0032 
REWIND 41 
0033 
REWIND 31 
0034 
REWIND 32 
0035 
REWIND 22 
0036 
REWIND 42 
0037 
REWIND 25 
0038 


B 
8lANK=606060606060 
0039 


8 
AST=545454545454 
0040 


8 
PLUS=202020202020 
0041 


8 
PER 100=333333333333 
0042 
0(1)=0.0 
0043 
TII1EW=I.0E-I0 
0044 
0045 
40 
FORI1ATC17HCONTROL DATA RUN 914.6X,EtO.3.tlAl) 
0046 


41 
FORI1AT(14HCONSTANTS RUN 12.7H CYCLE 14.1X.IP3E12.5.0P4F4.2) 
0047 


42 
FORI1ATCI7HDATA CARDS CYCLE 14.10X.I0HRUN NUMBERI3.10X.5HSET AI 
0048 


42 
lCI4.1P5EI2.5.1PElt.4.I2.I3)) 
0049 


43 
FORI1AT(17HDATA CARDS CYCLE 14. lOX. 10HRUN NUMBERI3.10X.5HSET 81 
0050 


43 
1C I 4. 1 P5E 1 2 • 5. 1 PE 1 1 • 4. I 2. I 3 ) ) 
0 05 I 
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44 
FORMAT(lH .SH MOOE=14/SH 
NMC=14/SHNZONE=14/SH 
NQS=!4/SH 
MCL=!41 
0052 
44 
ISH 
RS=IPEll.4/SH TIME=Elt.4/SH 
DT=Ell.4/SH 
CS=Ell.4/SH 
CR= 
0053 
44 
2Ell.4/SH 
BR=El1.4/7H TIMEL=lPEll.4) 
0054 
4S 
FORMAT(lHl.214.IPEI2.3.15.1PSEI2.3/(14.1P12ElO.3) 
0055 
47 
FORMAT(IH .19HZONE SPLIT AT NMC=14.5X.5HTIME=IPE11.4.SX.5HZONE=14. 
0056 
47 
119) 
0057 
48 
FORMAT(IHl.30X.4HRUN=14.15X.20HLOG POWER (ERGS/SEC)/6H CYCLEoX.4HT 
0058 
48 
2IME.3X.12.13X.12.13X.12.13X.12.13X.12.13X.12.X.11HTOTAL POWER.3X. 
0059 
48 
36HRAOIUS/15X.76AI) 
0060 
49 
FORMAT(lH .1SHREZONED AT NMC=I4.5X.SHTIME=IPE11.4.5X.5HZONE=14. 17) 
0061 
50 
FORMAT(14.IPSEI3.5) 
0062 
NSZ=l 
0063 


REWIND 15 
0064 


REWIND 16 
0065 


READ INPUT TAPE 5.40.NR.MODE.NMC.NLP.NDD.LZ.NZONE.NCS.MCL.TIMEL 
0066 


1.8LANK.AST.TEE.PLUS.PERIOD.OASH.ECUAL.PINUS.FFF.UUU.PPP 
0067 


PRINT 9321. MCL.TIMEL 
0068 


9321 
FORHAT(34H NORMAL TERMINATION CONDITION NMC=I4. IOX.5HTIME=IPEI2.3) 
0069 


WRITE OUTPUT TAPE 15.12345.NR 
0070 


12345 FORMAT(lHIIIIIIIIIII11H .40X. 
32HR W HILLENDAHL PALO ALTO 201 
0071 
123451111H .40X.33HPROOUCTION OUTPUT LIST RUN NUMBERI4) 
0072 


READ INPUT TAPE 5.41.NR.NMC.TIME.OT.FLEX.CS.CR.BR.RS 
0073 
NSTART=NHC 
0074 
DTMIN=DT 
0075 
NEG=NMC 
0076 
MCOC=NMC+50 
0077 
NMA6=HCDC 
0078 
JXC=O 
0079 
NHCS=NHC 
0080 
TIMES=TIME 
0081 


FLIX=FLEX 
0082 
NSTOP=O 
0083 
NRZ=50 
0084 


NB=14 
0085 


IPS=O 
0086 
LZHI=LZ-I 
0087 
LZH2=LZ-2 
0088 
IF(MOOE-I)I. I. 1000 
0089 


CALL ENTRY 
0090 


WRITE OUTPUT TAPE 6.44 
.MODE.NMC.NZONE.NQS.MCL.RS.TIME.DT.CS.CR. 
0091 
IBR.TIMEL 
0092 


FBR=R(100) 
0093 
RTEST=RCI) 
0094 


NPL=NHC 
0095 
GO TO 1036 
0096 
1000 
READ INPUT TAPE 5.42.NMC.NR.(NHC.RCK).U(K).V(K).QeKl.TCKl.P(K).NR. 
0097 
1000 lK.K=I.100) 
0098 
REAO INPUT TAPE 5.43.NMC.NR.eNMC.RzeKl.R2eKl.ZHAseKl.DReK).RHozeKl 
0099 


1.TReK).NR,K.K=I.IOO) 
0100 


READ INPUT TAPE 5.50.NR.TIME.TEE.WKT.TDN2,TD02 
0101 


READ INPUT TAPE 5.50.NR.SCALE.EQUAL.Q7I.Q47.TIMEW 
0102 
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~B=NHC+l 
0103 
LZ=100 
0104 
NLP=~LP+25 
0105 
DTHI~=DT.O.1 
0106 
CALL REJUST 
0107 
(PS=O 
0108 
~Pl=~HC 
0109 
1034 WRITE OUTPUT TAPE 6.1035.DT.DR(1).CK.R(K).RHOZ(K).V(K).TCK).U(K).Q 
0110 
1(K).ZHASCK).DRCK).K=I.100) 
0111 
1035 FORHATCIHI.IOHl~PUT OATA.IOX.3HOT=EI2.4.5X.3HOR=E12.4/3H 
K.2X.4HR 
0112 
ltL).8X.7HRHOZCK1.5X.4HVCK).8X.4HTCK).8X.4HUCL).8X.4HQCKl.8X.7HZMAS 
0113 
ICK).5X.5HORCK)/(14.1P8EI2.4)) 
0114 
0115 


WRITE OUTPUT TAPE 6.44 
.MOOE.~MC.~ZO~E.NQS.MCL.RS.TIME.OT.CS.CR. 
0116 


IBR.TIMEL 
0117 


FBR=RCIOO) 
0118 


RTEST=R(l) 
0119 


1036 
NTAPE=15 
0120 
GO TO 3000 
0121 


CHASTER CYCLE RE-ENTRY POINT 
0122 


2000 TIHE=TIME+OT 
0123 
CROUT INE TO CHANGE NUMBER OF ZONES IN USE .... LIMIT 100 ZONES 
0124 
10 
DO 14 
J=I.l00 
0125 
K=IOI-J 
0126 


IFCUCK)-I.OE+02) 11.12.12 
0127 


11 
IFCTCK)-2.0-TCI00)) 14.14.12 
0128 
12 
LZ=K+6 
0129 


13 
GO TO 16 
0130 


14 
CONTINUE 
0131 


15 
LZ=100 
0132 
16 
IF (L Z -I 00) 18. 18. 17 
0133 
17 
LZ=100 
0134 
18 
LZMJ=LZ-I 
0135 
19 
LZM2=LZ-2 
0136 


20 
LZPl=LZ+l 
0137 
21 
LZP2=LZ+2 
0138 
22 
FLEX=FLEX+FLOX*OT 
0139 
23 
IF(~MC) 33.24.26 
0140 
24 
NTAPE=15 
0141 
CNCLZ)=O.O 
0142 
TFLUX=FLOX 
0143 
FLEX=O.O 
0144 


25 
GO TO 732 
0145 


26 
IF(~MC- 3) 
27.27.29 
0146 


27 
NTAPE=6 
0147 


28 
GO TO 32 
0148 


29 
IF(TIME-TIMEV) 
33.30.30 
0149 
30 
TIMEW=TIMEV.C10.0 •• (1.0/18.0)) 
0150 
IFCNHC- 4) 
31.29.31 
0151 


31 
NTAPE=15 
0152 
LZP2=Nt 
0153 
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I 
32 
CN(LZ1=(TIME-TIMES1/FLOATF(NMC-~Sl 
0154 
TFLUX=(FLEX-FLIX1/(TIME-TIMESl 
0155 


732 
CALL CGSPO 
0156 


NOP-NaP 
0157 


FLIX=FLEX 
0158 


NMCS=NHC 
0159 
TIMES=TIME 
0160 


LZP2=LZ+2 
0161 
33 
CONTINUE 
0162 
0163 


494 
IF(NHC-HCOC) 
414.409.409 
0164 
409 
HCOC=HCOC+50 
0165 
RE~INO 16 
0166 
HOOE=HOOE+I 
0167 


~RITE OUTPUT TAPE 16.40.NR.HOOE.NHC.NLP.NEG.LZ.NZONE.NQS.HCL. 
0168 
ITIHEL.BLANK.AST.TEE.PLUS.PERIOO.OASH.EQUAL.PINUS.FFF.UUU.PPP 
0169 


~RITE OUTPUT TAPE 16.41.NR.NHC.TIHE.OT.FLEX.CS.CR.BR.RS 
0170 


~RITE OUTPUT TAPE 16.42.NHC.NR.CNHC.R(K).UCK).VCK).QCK).TCK). 
0171 


1 PC K ) • NR. K. K = 1 • 1 00) 
0 I 72 


WRITE OUTPUT TAPE 16.43.NHC.NR.CNHC.RZ(K).R2CK).ZHASCK).ORCK). 
0173 
lRHOZCK).TRCK1.NR.K.K=I.100) 
0174 
WRITE OUTPUT TAPE 
16.50.NR.TIHE.TEE.WKT.TON2.T002 
0175 


~RITE OUTPUT TAPE 16.S0.NR.SCALE.EQUAL.Q71.Q47.TIHEW 
0176 
HOOE=HOOE-I 
0177 
0178 
C REZONE 
S~ITCH 
0179 
414 
RTEST=RCLZ-Sl 
0180 
NZTS=LZ-l1 
0181 
IFCHCP-3) 
415.415.4150 
0182 
C 
4150 
ALL OPTICALLY THICK CASES 
0183 
4150 
NZTS=HCP-4 
0184 


PHT=2.5S 
0185 


4153 
IFCNZTS-9) 
427.417.417 
0186 
C 
415 
ALL TRANSPARENT CASES 
0187 


415 
IFCT(4)-2.0.TON2) 
4170.4170.416 
0188 
C 
416 
HIGH ALTITUDE EARLY PHASE 
0189 


4 16 
00 4 1 61 
L = 9. L Z 
0 I 90 
IFCTCL)-I.S.TON2) 
4160.4161.4161 
0191 
4160 
NZTS=L-l 
0192 


PHT=I.S 
0193 
IFCNZTS-9) 
427.417.417 
0194 


4161 
CONTINUE 
0195 
C 
4170 
LATE TIME TRANSPARENT CASE 
0196 
4170 
00 4140 
J=11.LZ 
0197 
IFCRCJ)-WCl.9)) 
4140.4141.4141 
0198 


4141 
NZTS=J-5 
0199 


PHT=2.8 
0200 
IFCNZTS-9) 
427.417.417 
0201 
4140 
CONTINUE 
0202 
417 
00 426 
K=9.NZTS 
0203 
4013 
IFCNHC-NRZ) 
427.427.419 
0204 
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419 
422 
423 
424 
425 
9425 


IFCCRCKl-RCK-2ll-0.0a.RTESTl 
422.422.426 


IFCABSFCTCKl-TCK-lll-0.30.TCK» 
423.423.426 


IFCABSFeYCKl-YCK-l»-O.7-YCKl) 
424.424.426 
IFCZHASCK)+ZHASCK-l)-PHT.ZHASCK-2l) 
425.425.426 


IF(QCK)-0.05.PCK» 
9425.9425.426 


LZP2-K 
NRZ=NHC 
JXC=JXC+l 
WRITE OUTPUT TAPE22.49.NMC.TIME.LZP2 .JXC 
WRITE OUTPUT TAPE 6.49.NHC.TIHE.LZP2 .JXC 
CAll REZONE 
lZ=LZ-l 
lZH2=LZ-2 
LZH1=lZ-t 
lZPI=LZ+l 
lZR=LZR-2 
llP2=lZ+2 
NRZ=NHC+2 
GO TO 427 
426 
CONTINUE 
C 
OPTICAL SPLIT TEST 
427 
NS=LZ-5 


428 
429 


430 
431 
432 
433 
434 
435 


436 


437 


440 
441 
442 
443 


444 


KS=O 
IFCNMSP-NMC) 
428.428.448 


00 440 
J=3.NS 


K=LZ-J 
IFC(RCK)-RCK-l))-0.08.RCLZ)) 
430.430.442 
IF(TAUCK.l)-1.7) 
440.440.431 
KS=KS+i 
IFCKS-I) 
448.433.448 


-JFCQ(Kl-PCK» 
448.434.434 
IF(UCK1-l.OE+05) 
448.435.435 


KST=K 
KSP=K+5 
00 436 
J=KST.KSP 


IF((RCJ1-RCJ-l11-200.0) 
436.436.437 


CONTINUE 
GO TO 448 
KZ2=J 
IF(TCJ)-6.0E+041 
443.448.448 
CONTINUE 
GO TO 448 
KZ2=K 
NHSP=NMC+3 
JXO=JXO+l 
CALL SPLIT 
WRITE OUTPUT TAPE 
6.47.NMC.TIME.KZ2.JXO 
WRITE OUTPUT TAPE 22.47.NMC.TIME.KZ2.JXO 
LZ=LZ+l 
LZf11=LZ-l 
LZPl =LZ+ 1 
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0205 
0206 
0207 
0208 
0209 
0210 
0211 
0212 
0213 
0214 
0215 
0216 
0217 
0218 
0219 
0220 
0221 
0222 
0223 
0224 
0225 
0226 
0227 
0228 
0229 
0230 
0231 
0232 
0233 
0234 
0235 
0236 
0237 
0238 
0239 
0240 
0241 
0242 
0243 
0244 
0245 
0246 
0247 
0248 
0249 
0250 
0251 
0252 
0253 
0254 
0255 


LZP2-LZ+2 
0256 
L~2.LZ-2 
0257 
0258 
448 
IF(NHC-NSTART) 
449.451.449 
0259 
449 
IF(~-NSTA.T-l) 450.451.450 
0260 
450 
IF(~-NOO) 453.451.453 
0261 


451 
KZ2=-451 
0262 
452 
CALL OIAGNS 
0263 
453 
NMC=N~+1 
0264 
0265 
C 
0266 
IF(NSTOP) 
35.35.36 
0267 


35 
IFnIl'1E-6.0) 34.52.52 
0268 


34 
IF[SENSE SWITCH 1) 
52.65 
0269 


52 
"ooE="OOE+l 
0270 
NSTOP= 1 
0271 
~L="CL+3000 
0272 
CALL OIAGNS 
0273 
GO TO 28 
0274 


36 
PUNCH 40.NR.HODE.NHC.NLP.NEG.LZ.NZONE.NQS.MCL.TIHEL 
0275 
1.BLANK.AST.TEE.PLUS.PERIOD.DASH.EQUAL.PINUS.FFF.UUU.PPP 
0276 
PUNCH 41.NR.NHC.TIME.DT.FLEX.CS.CR.8R.RS 
0277 
PUNCH 42. NHC.NR.CNHC.RCK).UCK).VCK).QCK).TCK).P(K).NR.K.K=I.100) 
0278 
PUNCH 43.NMC.NR.CNMC.RZ(K).R2(K).ZHASCK).DRCK).RHOZCKl 
0279 
I.TRCK1.NR.K.K=I.100) 
0280 
PUNCH 50.NR.TIHE.TEE.WKT.TDN2.TD02 
0281 
PUNCH 50.NR.SCALE.EQUAL.Q71.Q47.TIMEW 
0282 
END FILE 15 
0283 


END FILE 25 
0284 
REWIND 25 
0285 
CALL COpy C 25. 15) 
0286 
END FILE 15 
0287 


END FILE 22 
0288 
REWIND 22 
0289 
CALL COPY(22.15) 
0290 


END FILE IS 
, 
0291 
WRITE OUTPUT TAPE 6.37.JXC 
0292 


37 
FORHAT(IH .13HREZONE CALLEDI4.6H TIMES) 
0293 
REWIND 16 
0294 


NR=9999 
0295 
WRITE OUTPUT TAPE 41. 1499.NR.KZ6 
0296 
wRITE TAPE 31.NR.KZ6 
0297 


1499 FORMAT(214) 
0298 
WRITE OUTPUT TAPE 32. 45.NR.NHC.TIHE.LZ.P(99).TC99).RHOZC99).POWER 
0299 
2.P47.FBR.CK.RCK).UCK).PCK).URCK).TCK).FOSCK).FISCK).DHCK).QCK). 
0300 
3FZ(K).ECK).HPPPCK).K=I.LZ) 
0301 
IFC....c:-NHAS) 
1502.1502.1501 
0302 
1501 
READ INPUT TAPE 
16.40.NR.HOOE.NHC.NLP.NEG.LZ.NZONE.NQS.HCL. 
0303 
ITIHEL.BLANK.A$T.TEE.PLUS.PERIOD.OASH.EQUAL.PINUS.FFF.UUU.PPP 
0304 
READ INPUT TAPE 
16.41.NR.NHC.TIHE.OT.FLEX.CS.CR.BR.RS 
0305 
READ INPUT TAPE 
16.42.NHC.NR.CNHC.RCK).U(K).V(K).QCK).TCK). 
0306 
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IPCK).NR.K.K=l.lOO) 
0307 


READ INPUT TAPE 
16.43.NHC.NR.(NHC.RZCK).R2CK).ZHASCK).DRCK). 
0308 
lRHOZCK).TRCK).NR.K.K=1.100) 
0309 
READ INPUT TAPEI6.50.NR.TIHE.TEE.WKT.TDN2.TD02 
0310 
READ INPUT TAPEI6.S0.NR.SCALE.EQUAL.Q71.Q47.TIMEW 
0311 
PUNCH 40.NR.HODE.NHC.NLP.NEG.LZ.NZONE.NQS.HCL.TIMEL 
0312 
1.BlANK.AST.TEE.PLUS.PERIOD.DASH.EQUAL.HINUS.FFF.UUU.PPP 
0313 
PUNCH 41.NR.NHC.TIHE.DT.FLEX.CS.CR.BR.RS 
0314 
PUNCH 42. NHC.NR.CNHC.RCK).UCK).VCK).QCK).TCK).PCK).NR.K.K=I. 100) 
0315 
PUNCH 43.NHC.NR.(NHC.RZeK).R2eK).ZHASCK).DRCK).RHOZCK) 
0316 
1.TRCK).NR.K.K=I.100) 
0317 
PUNCH 50.NR.TIHE.TEE.WKT.TDN2.TD02 
0318 
PUNCH SO.NR.SCALE.EQUAL.Q71.Q47.TIHEW 
0319 


END FILE 16 
0320 
1502 
CALL UNLOADCI6) 
0321 
CEASE=1.0E+29 
0322 


END FILE 32 
0323 


END FILE 42 
0324 
REWIND 42 
0325 
CALL COPY(42.32) 
0326 
CALL UNlOAD(42) 
0327 
0328 
C 
EXTRA COpy OF USER TAPE 32 INFO COPIED AS FILE 2 OF TAPE 15 
0329 


END FILE 32 
0330 
REWIND 32 
0331 
CALL COPYC32.15) 
0332 


END FILE 15 
0333 
CALL COpy (32.1S) 
0334 
EN() FILE 15 
0335 
CALL UNLOAD(15) 
0336 
CALL UNLOAD(32) 
0337 


END FILE 41 
0338 
CALL UNlOAD(41) 
0339 


END FILE 31 
0340 
CALL UNLOAD(31) 
0341 
IS00 
CALL EXIT 
0342 


65 NT[=O 
0343 
IFCNMC-HCLl 
64.64.52 
0344 
64 
IF(TIHE-TIHEL) 
66.66.52 
0345 
C SET RETARDED VARIABLES 
0346 


66 
DO 71 
K=I.I00 
0347 


67 
VRCK)=V(K) 
0348 


68 
ERCK)=ECK) 
0349 


69 
PRCK)=P(K) 
0350 
70 
VCK.4l=T(K) 
03S1 
WCK.6)=VCK) 
0352 
W(K.7,=PCK) 
0353 
WCK.8)=QCK) 
0354 
WCK.9)=ECK) 
0355 


71 
DIVFRCK)=DIVFACK) 
0356 
KRT=O 
0357 
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72 
IF(DT-50.0-DTHINl 
74.74.73 
73 
DTHIN=DT/50.0 
74 
CALL HYDRO 
KZ4=KZ4 
IFCU(99l-1.0E+03) 290.290.52 
C 
INITIAL TEMPERATURE EXTRAPOLATION 
290 
CONTINUE 
C 
ITERATION CYCLE RE-ENTRY POINT 
C DERIVATIVES 
WITH RESPECT TO TEMP 
OP 
3000 
00 301 
K=I.lZ 
3001 
IFCTCK» 
3002.3002.300 
3002 
KZ2=-3000 
CAll OIAGNS 
GO TO 52 


300 
TP(K)=1.02-TCK) 
301 
TSCK)=TP(K) 


303 
CAll STATE 


304 
00 310 
K=I.lZ 


305 
PPCK)=PCK) 
306 
EPCK)=ECK) 
309 
TMCK)=0.98-TCK) 
310 
TSCK)=TMCK) 


311 
CALL STATE 
312 
00 
318 
K=I.LZ 


313 
PMCK)=PCK) 


314 
EMCK)=ECK) 


315 
DPTCK)=(PPCK)-PMCK»/CTPCK)-TMCK» 


316 
DETCK)~CEPCK)-EMCK»/CTPCK)-TMCK» 


318 
TSCK)=T(K) 
320 
CALL STATE 


C CALL RADIATIvE PROPERTIES ROUTINE 
CALL SWABZ 
501 
LZR=LZR 
MCP=MCP 
KZ1=KZl 
KZ2=KZ2 
KZ3=KZ3 
KZ8=KZ8 
KZ9=KZ9 


321 
NCW=NCW 


322 
DTAUCLZP1.1)=0.0 


C 
INTENSITY INTEGRATION 
500 
CAll FLUXS 


502 
IFCIPS) 
503.503.530 
503 
IPS=IPS+l 


504 
TIHE=TIME-DT 
505 
GO TO 2000 
C 
TEMPERATURE TEST BYPASS ON FIRST GUESS 


530 
IFCNTI) 
540.540.532 
C TEMPERATURE TEST 
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DE 


0358 
0359 
0360 
0361 
0362 
0363 
0364 
0365 
0366 
0367 
0368 
0369 
0370 
0371 
0372 
0373 
0374 
0375 
0376 
0377 
0378 
0379 
0380 
0381 
0382 
0383 
0384 
0385 
0386 
0387 
0388 
0389 
0390 
0391 
0392 
0393 
0394 
0395 
0396 
0397 
0398 
0399 
0400 
0401 
0402 
0403 
0404 
0405 
0406 
0407 
0408 


532 
DO 536 K=l.LZ 
0409 


533 
IF(ABSF(OTt1(K1IT(Kll-0.1 1 .536.536.540 
0410 
536 
CONTINUE 
0411 
C t1AIN CYCLE COt1PLETEO -- 
RETURN TO 2000 
0412 


531 
GO TO 2000 
0413 
540 
NTI=NTJ+I 
0414 
C ENERGY EQUATION 
BLOCK 600 
0415 
600 
IF(NTB-3) 601.601.2000 
0416 
601 
DO 617 K=t.LZ 
0417 
602 
L=K 
0418 
603 
N=llO+K 
0419 


613 
OECK1=CECK1-ER(K»/DT 
0420 
614 
PACK1=CPRCK)+PCK»/2.0 
0421 


615 
RESDUECKl=DECK)+CPACK)+QCK».CVCK1-VRCK)/DT 
0422 


625 
1+(DIVFACK)+DIVFRCK»/2.0 
0423 
617 
CONTINUE 
0424 
C 
TE~RATURE ITERATION OF ENERGY EQUATION 
0425 
700 
IFCLZR-4) 
703.703.701 
0426 
701 
CALL COEFF 
0427 
702 
GO TO 829 
0428 
C 
NO ITERATION OF RADIATION 
0429 
703 
00 704 K= 1. LZ 
0430 
704 
DTHCK1=-RESDUECK)/C(OETCK)/DT)+CDPTCK)/C2.0.DT».(V(K)-VRCK»)) 
0431 
705 
IFCNHC-NlP) 
906.706.906 
0432 
706 
NLP=NLP+50 
0433 


PRINT 707.NMC.TIME.DT.LZ 
0434 
707 
FORHATC23H PROGRESS REPORT 
NHC=14. 5X.5HTIME=IPEI2.3.5X.I1HTIME 
0435 
707 
I STEP=lPEI2.3.10X.5HNORAD(4111) 
0436 
708 
GO TO 906 
0437 
0438 
C 
RADIATIVE ITERATION 
0439 


829 
CONTINUE 
0440 
830 
IF(NMC-NlP) 
987.831.987 
0441 


831 
NLP=NLP+50 
0442 
PRINT 832.NMC.TIME.DT.LZ 
0443 


832 
FORt1AT(23H PROGRESS REPORT 
NMC=14. 5X.5HTIME=IPE12.3.5X.l1HTIME 
0444 


832 
I STEP=IPEI2.3.10X.8HRADHYDROI4111) 
0445 
C 
0446 
C 
COEFFICIENTS NOW KNOWN-SOLUTION FOR DTMeK) 
0447 


987 
LZJ=lZR+l 
0448 
00988 
K=LZJ.LZ 
0449 
988 
DTI1( K) =-RESDUECK )/( (OETCK) lOT) +( DPTCK) I( 2. O.DTl).( V( K) -VRCK») 
0450 
900 
OTH(LZR)=(ENCLZR)-ENCLZR-l»/CCNCLZR)-CN(LZR-I» 
0451 
902 
DTH(LZR-l)=(EN(LZR-I)-CNCLZR-I).DTMCLZR» 
0452 
DO 
905 
J=3.LZR 
0453 
K=LZJ-J 
0454 
905 
OTt1CK'=ENCK1-ONCK).OTt1CK+2)-CNCK).OTt1CK+t) 
0455 
906 
CONTINUE 
0456 
C STABILITY CHECK 
0457 
907 
00 930 
K=I.LZ 
0458 
908 
IF(OTI1(K») 
909.930.910 
0459 
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909 
IFCABSFCDTHCK»-0.5.TCK») 
930.930.918 
0460 
910 
IFCTCK)-2.0E+OS) 
912.911.911 
0461 


911 
IFCABSF(DTM(K»-0.9.T(K» 
~30.930.918 
0462 


912 
IFCABSF(DTHCK»-S.O.TCK» 
930.930.918 
0463 


918 
KZ2=-920 
0464 


919 
KZ5=K 
0465 
920 
CALL DIAGNS 
0466 


921 
KRT=KRT+l 
0467 


922 
IFCKRT-3) 
923.933.933 
0468 


923 
00 927 
L=I.LZ 
0469 
U(Ll=VCL.3) 
0470 
T(L)=VCL.4) 
0471 
RCL)=VCL.S) 
0472 
VCL)=VCL.6) 
0473 
PCL)=VCL.7) 
0474 
QCL)=VCL.8) 
0475 
ECL)=VCL.9) 
0476 
R2CL)=RCL) •• 2 
0477 


927 
CONTINUE 
0478 


928 
VRITEOUTPUT TAPE 6.929.NMC.KRT.K 
0479 


PRINT 929.NMC.KRT.K 
0480 


929 
FORMAT(4I8) 
0481 
NTI=O 
0482 
GO TO 74 
0483 
930 
CONTINUE 
0484 


C MODIFY TEMPERATURE 
0485 


933 
00 934 
K=I.LZ 
0486 


934 
HK)=TCK)+DTMCK) 
0487 


935 
GO TO 3000 
0488 


936 
END 
0489 
• 
FORTRAN 
0490 
CF1SVZ 
0491 
SUBROUTINE SWABZ 
0492 


C 
STANDARD DIMENSION AND COMMON STATEMENTS 
AS IN MAIN PROGRAM 
0493 
D(1)=0.0 
0494 
0495 
DO 126 
K=I.LZ 
0496 
00 135 
N=I.NSZ 
0497 
TSCK)=TSCK)/11606.5 
0498 
ETACK)=1.293E-03.V(K) 
0499 
ALPHA=9.0+EXPFC-TIME/I.OE-03.CETACK) •• 2» 
0500 
BETA=CC6.0E-18.VCK» •• O.5).CTSCK) •• 0. I) 
0501 
300 
ZICK.I)=EXPFC-I.S.CRCK)-RCK-I» 
ICCETACK) •• I.S)/CCTSCK) •• ALPHA) 
0502 
300 
2+BETA) 
+0.2.CETACK) •• 1.91)/CTSCK) •• 2.73) 
0503 
300 
3+0.023.CETACK) •• 1.8).CTSCK) •• 0.2S)+I.OE-07.CETACK) •• 2.0)1 
OS04 
300 
4CCTSCK) •• C-S.O»+4.0E-14.VCK»» 
0505 
ZICK.I)=ZICK.I).0.99999998 
OS06 
ZICK.2)=ZICK.I) 
0507 
HHFP=CETACK) •• I.SJ/CCTSCK) •• ALPHA)+BETA) 
0508 
1+0.2.CETACK) •• 1.91)/CTSCK) •• 2.73) 
OS09 
2+0.023.CETACK) •• 1.8).CTSCK) •• 0.2S) 
0510 


138 


3+1.0E-07·(ETA(~) •• 2.0).(TS(~) •• 5.00) 


HHFF·9.3E-06.(ETA(~) •• 2.0).(TS(~)·.3.5) 


IF(TS(~)-IO.O) 301.298.298 
298 
IF(HHFF-HMFP) 
299.301.301 


299 
t-I1FP -Ht1F F 


301 
OTAU(~.I)=1.5.(R(~)-R(~-1»/HMFP 
302 . A(~.N)=EXPFC-OTAUC~.N» 
312 
8CC~.N)=5.67E-OS.CTC~) •• 41 


TSC~)=TSC~).11606.5 


135 
CONTINUE 
126 
CONTINUE 


C 
SOURCE FUNCTION INTERPOLATION 
DO 201 
~=l.LZ 
DO 202 
N=l.NSZ 
311 
8B(K.N)=CC((CRC~+I)-RC~l).T C~) •• 4)+CCRCK)-RCK-l».T CK+l) •• 4»1 


311 
1(RCK+l1-RCK-l»».5.67E-05 
BBCK.2)=BBCK.l) 
506 
SCK.Nl=CBCCK.Nl-BCCK+l.N»/CDTAUCK.N)+OTAUCK+l.N» 
507 
WCK.Nl=I.O-ACK.N)-ACK.Nl.OTAUCK.N) 
508 
IF(W(K.Nl-l.0E-04) 
509.202.202 
509 
WCK.Nl=0.5.COTAU(K.N) •• 2) 


202 
CONTINUE 
201 
CONTINUE 


C 
SET OPTICAL INDICES 
Kf1=O 
501 
KD=O 
502 
KN=O 
503 
TAUCLZ.l)=DTAU(LZ.l) 


511 
DO 
522 
J=I.LZf11 


512 
K=lZ-J 


513 
TAUCK.ll=TAUCK+l.1)+OTAUCK.l) 


514 
IF(TAUCK.l1-5.0E-06) 
522.515.515 
515 
KN=KN+l 


516 
IFCKN-l) 
518.517.518 
517 
lZR=K+3 


518 
IFCTAU(K.l1-.44) 522.519.519 
519 
KD=KO+l 
520 
IFCKD-l) 
622.523.622 
523 
FBR=RCK) 


622 
IFCTAUCK.l)-1.0) 
522.623.623 
623 
KM=KM+l 
624 
IFCKf1-1) 
522.625.522 
625 
MCP=K 


522 
CONTINUE 
525 
IFCLZR-LZ+1) 530.530.529 


529 
LZR=lZ-1 
530 
CONT INUE 
IF(KO) 
521.521.524 
521 
FBR=O.O 
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0511 
0512 
0513 
0514 


051~ 
0516 
0517 
0518 
0519 
0520 
0521 
0522 
0523 
0524 
0525 
0526 
0527 
0528 
0529 
0530 
0531 
0532 
0533 
0534 
0535 
0536 
0537 
0538 
0539 
0540 
0541 
0542 
0543 
0544 
0545 
0546 
0547 
0548 
0549 
0550 
0551 
0552 
0553 
0554 
·0555 
0556 
0557 
0558 
0559 
0560 
·0561 


'1(P-O 
0562 
524 
CONTINUE 
0563 
C 
SET CUTOFF ~AVELENGTH 
C~L(K) 
0564 
CWHAX=O.O 
0565 


KC~HAX=O 
0566 
C~L(HCP.I)=O.O 
0567 


NC~=HCP 
0568 
TR01CP)=TCH(P) 
0569 
0570 
540 
00 
565 
K=I.LZ 
0571 


541 
IF(T(K)-TDN2) 
546.542.542 
0572 
0573 
C 
OPACITY DUE TO ATOMIC SPECIES 
0574 


542 
CWL(K.l)= 700.0-EXPFC-0.36.CTCK)-11606.5)/11606.5) 
0575 


8544 
DTH(K)=212121212121 
0576 


IFCC~LCK.l)-275.0) 
545.560.560 
0577 


545 
C~L(K. 1 )=275. 0 
0578 
GO TO 560 
0579 
546 
IF(TCK)-TD02J551.547.547 
0580 
0581 
C 
OPACITY DUE TO N2 HOLECULE 
0582 


8 547 DTHCK)=454545454545 
0583 
TCN2=5.0E-07-1.293E-03-VCK) 
0584 


549 
CWLCK.l)=1140.0.C(CRCK)-RCK-I»/CYCKJ-l.293E-03» •• 0.1I) 
0585 
l-CI.0-EXPFC-TIME/TCN2» 
0586 
IFCCWLCK.l)-IOOO.O) 
550.560.560 
0587 
550 
C~LCK.l)=IOOO.O 
0588 
GO TO 560 
0589 
0590 
C 
OPACITY DUE TO 02 MOLECULE 
0591 


551 
TC=4250.0-271.0.LOGFCI.293~-03. YCK» 
0592 
TC02=3.0E-07.1.293E-03-YCK) 
0593 


8 552 DTHCK)=676767676767 
0594 


553 
IFCTCK1-TC1554.554.556 
0595 


554 
CWLCK.l)=1500.0+TSCK)-CO.163+0.0745.LOGFCDRCK).RHOZCK)11.293E-03» 
0596 


554 
l-C].0-EXPFC-TIME/TC02J) 
0597 


555 
GO TO 557 
0598 


556 
TS8=TSCK).C.0647-LOGFCRCKJ-RCK-IJ)-0.25-.109.LOGFC1.293E-03.V(K») 
0599 
CWL(K.l1=3500.0+TSB-C2000.0+TS8)-EXPFC-TIME/TC02) 
0600 


557 
IFCCWLCK.])-1500.0) 558.560.560 
0601 


558 
C~L(K.I )=1500. 0 
0602 
560 
IFCK-MCPl 
565.565.561 
0603 


561 
IFCCWLCK.ll-CWMAX) 
565.565.563 
0604 


563 
CWHAX=CWLCK.l) 
0605 


564 
KC~"AX=K 
0606 


565 
CONTINUE 
0607 


566 
NCW=KC~HAX 
0608 
0609 
IFCLZR-NC~) 
603.603.604 
0610 
603 
LZR=NCW+ 1 
061 I 
C 
CALCULATE ZO VALUES 
0612 
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604 
ZO(I.I)-ZI(I.I) 
KO=O 
605 
DO 670 
L=2.LZ 
TRCL )=TCL) 
IF(TCLl-9.0E+041 
607.669.669 


607 
KO=KO+I 
IF(KO-I) 
609.609.609 
609 
P1S=L-l 
609 
IF(DTAUCL.ll-l.O) 
610.669.668 


C 
SELECT TRCll= TEMPERATURE OF RADIATION TRAVERSING ZONE L 


610 
LIP1=l-1 


611 
T AUSUM=O. 0 


615 
00 630 NN=I.LIM 


616 
P1=l-NN 


617 
TAUSUP1=TAUSUM+DTAUCM.I) 
620 
IF(TAUSUM-0.7) 
630.627.627 


627 
TRCL) =T(M) 


628 
IFCP1-2), 630.630.638 
630 
CONTINUE 
IFCT(Ll-2500.0) 
631.632.632 
631 
TRCL)=TRCL-l) 
GO TO 637 
632 
EPTSUM=O.O 
EPSUM=O.O 
00 635 
K=3. LlP1 
EPSUP1=EPSUP1+DTAUCK.l) 
EPTSUP1=EPTSUM+TCK)-OTAUCK.I) 
IFCT(K)~2500.0) 
636.636.635 
635 
CONTINUE 


636 
TRCL)=EPTSUM/EPSUM 


637 
"'=2 


638 
I Z=P1+ 1 
639 
DO 646 
1=IZ.l 
640 
IFCTRCL)-CWLC].I)-2.0E+08) 
641.641.644 
641 
ZICI.2)= EXPFCC-3.4IE-08-CWLCI.I)-TRC l )+2.5E-25-CCCWlCI.I)-TRCl 


641 
2» __ 3»_EXPFC-3.07E+13/CCCWLCI.I)-TRC L »--1.8») 


642 
GO TO 646 
644 
ZI(I.2)=7.61E+22/CCCWLCl.I)_TRC l »--3) 


646 
CONTINUE 
ZP1AX=Z J( 1'1+ 1.2) 
ZOCL.ll=ZlCH+l.2) 


647 
IF(L-P1-1) 669.663.650 
C 
CORRECTION FOR INTERP1EDlATE ZONES 
650 
IZ=f1+2 
651 
CWLTST=CWLCH+l.l) 


652 
00 660 
NZ=lZ.l 


653 
IFCCWlCNZ.I)-CWLTST) 
654.654.657 
654 
)F(NZ-L) 
660.655.668 
655 
ZOCL.I)=1.0 


656 
GO TO 660 


657 
ZOCL.l)=ZICNZ.2)/ZHAX 
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,0613 
0614 
0615 
0616 
0617 
0618 
0619 
0620 
0621 
0622 
0623 
0624 
0625 
0626 
0627 
0628 
0629 
,0630 
0631 
0632 
0633 
0634 
0635 
0636 
0637 
0638 
0639 
0640 
0641 
0642 
0643 
0644 
0645 
0646 
0647 
0648 
0649 
0650 
0651 
0652 
0653 
0654 
0655 
0656 
0657 
0658 
0659 
0660 
0661 
0662 
0663 


658 
CVLTST-CVL(NZ.l) 
0664 


659 
lHAX-ZICNZ.2) 
0665 
660 
CONTINUE 
0666 


663 
IF(Z~CL.l)-ZICL.I» 
670.668.668 
0667 


668 
Z~ ( L. 1 ) e Z I C L. 1 ) 
0668 
670 
CONT INUE 
0669 


671 
DO 676 
K=I.LZ 
0670 


672 
IF(OTAU(K.I)-I.OE-04) 
673.675.675 
0671 
673 
ACK. 1 )=OTAUCK. I) 
0672 


674 
GO TO 676 
0673 


675 
ACK. 1 )=1. O-ACK.I) 
0674 


676 
CONTINUE 
0675 


00 699 
K= I. LZ 
0676 
IFCSCK.l)-1.69E+38) 
767.767.677 
0677 
767 
HZ(K)=BBCK.I)-ACK+l.1) 
0678 
OZ(K)=2.0-SCK.I)-VCK+I.I) 
0679 
HPCK)=BCCK+l.1 )-ACK+l.l) 
0680 
IF(HZ(K)-OZ(K)-HP(K).O.OI) 
677.677.699 
0681 
677 
SCK.I)=O.O 
0682 


678 
BBCK.l)=BCCK+I.I) 
0683 


679 
BBCK. 2)=BCCK. I) 
0684 


699 
CONTINUE 
0685 
900 
IFCBCCI.I)-BC(2.1» 
915.916.916 
0686 


915 
sac I. 2)=BC( 1.1) 
0687 


916 
CONTINUE 
0688 
EMS= Z I C NCW .2) 
0689 


179 
RETURN 
0690 


END 
0691 
0692 
- 
FORTRAN 
0693 
CFIRJST 
0694 
SUBROUTINE REJUST 
0695 
C 
STANDARD DIMENSION AND COMMON STATEMENTS 
AS IN MAIN PROGRAM 
0696 
TD02=7100.0-EXPFCCO.43429-l0GFCRHOZ(99)/1.293E-03»/5.3) 
0697 
TDN2=15000.0-EXPF(CO.43429-LOGFCRHOZC99)/1.293E-03»/4.8) 
0698 


RETURN 
0699 


END 
0700 
- 
FORTRAN 
0701 
CFIRZNE 
0702 
SUSROU TINE RE ZONE 
07 03 
C 
STANDARD DIMENSION AND COMHON STATEMENTS 
AS IN MAIN PROGRAM 
0704 
0705 
H=LZP2 
0706 
Z,ETA= (T e M-l 1 +T(Ml )/(E (M-1 ) +E eM)) 
0707 
C 
ROUTINE TO COMBINE ZONES M AND M-I 
0708 


FKE= CZHASCH-2)-CUCM-2)--2)+ZMASCM-l1-((UCM-I)--2)+(U(M-2)--2) 
0709 
2)+ZHASCM)-CCU(M-l) •• 2)+(UCM) •• 2)+ZMASCM+I).(U(H) •• 2» 
0710 
FHV=CZHAS(H-2).U(M-2)+ZMASCM-I)-(UCM-l)+UCH-2»+ZHASCM)-CU(M-l) 
0711 
2+U(H»+ZMASCH+l).UCH» 
0712 


FKA=ZHAS(H-2)+ZHAS(H-I)+ZMASCH) 
0713 
FKB=ZHASCH-I)+ZHASCH)+ZHAS(H+l) 
0714 
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FKB-FKB/FKA 
FKE·FKE/FKA 
F"Y·F"Y/FKA 
FSR=CFKE.FKB1.CFKB+I.0)-FKB.CFMY·.2) 
IFCFSR) 
7700.7701.7701 
7700 
UC"-2)=UCM-2) 
UC"-tl=UCM) 
GO TO 7733 


7701 
FSR:,ICFSRIIO.5)/FKB 
UPLUS=(F"V+FSR)/CFKB+I.O) 
U"INS=CFMV-FSR1/CFKB+l.0) 
IFCU01-2)-UCM» 
7722.7711',7711 
7711 
UCM-2)=UPlUS 
UCM-l)=UMINS 
GO TO 7733 
7722 U(M-21=UMINS 
UCM-l1=UPLUS 


7733 
CONTINUE 


2 
ECM-ll=CZMASCM-I).ECM-I)+ZMASCM).ECM»/CZMAS(M-I)+ZMASCM» 


3 
TCM-l)=(ECM-I) 
'.ZETA) 
4 
ZMASCM-I)=ZMASCM-l)+ZMASCM) 
5 
0IVFA(M-l)=CR2CM)/ZMASCM-I».CFOSCM)-FISCM» 
5 
I-CR2(M-2)/ZMASCM-l».CFOSCM-2)-FISCM-2» 


7 
RCM-l)=RCM) 


B 
RzeH-ll=RzeM) 


9 
R2(M-I)=R2CM) 
10 
ORCM-ll=ORCM-l)+ORCMl 
II 
RRCM-t)=RRCM) 


13 
TRCH-I)=TRCM) 


14 
VCM-l)=CRCM-l) •• 3-RCM-2) •• 3)/C3.0.ZMASCM-I» 
15 
PCM-l)=CPCM-l)+PCM»/2.0 


16 
OCM-l)=COCM-l)+OCM»/2.0 
OTM01-1 )=0. 0 


C 
SHIFT IN EXTERIOR ZONES 
20 
00 50 
K=M.99 
21 
VCK)=VCK+l) 


22 
TCK)=TCK+l) 


23 
RHOZ(K)=RHOZCK+l) 


24 
UCK)=UCK+l) 
25 
PCK)=PCK+l) 


26 
RCK1=RCK+l) 


27 
R2eK)=R2CK+I) 


28 
RZ(Kl=RzeK+l) 


29 
OR(t('=ORCK+l1 
30 
O(K)=Q(K+I) 


31 
E(K)=E(K+l) 
32 
0IVFACK)=0IVFA(K+l1 


33 
URCK1=URCK+l) 
34 
ZHASeK)=ZHAS(K+l) 
35 
RR(K)=RRCK+I) 
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0715 
0716 
0717 
0718 
0719 
0720 
0721 
0722 
0723 
0724 
0725 
0726 
0727 
0728 
0729 
0730 
0731 
0732 
0733 
0734 
0735 
0736 
0737 
0738 
0739 
0740 
0741 
0742 
0743 
0744 
0745 
0746 
07<47 
0748 
0749 
0750 
0751 
0752 
0753 
0754 
0755 
0756 
0757 
0758 
0759 
0760 
0761 
0762 
0763 
0764 
0765 


DTHCK)=DTHCK+I) 
0766 


50 
TRCK).TRCK+l) 
0767 
0768 


C ADDITION OF NE~ ZONE 100 
0769 
51 
R(100)=R(99)+ORC99) 
0770 


52 
R2CIOO)=RCIOO) •• 2 
0771 
S3 
RZCIOO)=RCIOO) 
0772 
54 
RHOZCIOO)=RHOZC99) 
0773 


55 
RRCI001=1.0 
0774 
56 
DRCI001=DR(99) 
0775 


57 
ZMASCIOO)=CRHOZCIOO)/J.0).CCRZCIOO) •• 3)-CRZC99) •• 3») 
0776 


58 
UCIOO)=O.O 
0777 
59 
TCIOO)=Te99) 
0778 
60 
VeI001=V(99) 
0779 


62 
QCI001=0.0 
0780 
0781 
TCIOll=TeIOO) 
0782 


RETURN 
0783 


END 
0784 


• 
FORTRAN 
0785 
SUBROUTINE SPLIT 
0786 
C 
STANDARD DIMENSION AND COMMON STATEMENTS 
AS IN MAIN PROGRAM 
0787 
M=100-KZ2 
0788 
DO 30 
J=I.M 
0789 
K=IOI-J 
0790 
RCK)=RCK-l) 
0791 
DRCK)=DReK-l) 
0792 
RRCK)=RRCK-l) 
0793 


RHOzeK)=RHOZeK-l) 
0794 
RZCK)=RZeK-l) 
0795 
R2CK)=R2CK-l) 
0796 
UCK)=UCK-I) 
0797 
RRCK)=RReK-l) 
0798 
TRCK)=TReK-l) 
0799 
VCK)=VCK-l) 
0800 
PCK)=PCK-l) 
0801 
ZMASCK)=ZMASCK-l) 
0802 
QCK)=QCK-l) 
0803 
TCK)=TCK-I) 
0804 
ECK)=ECK-I) 
0805 
DIVFACK)=DIVFACK-l) 
0806 
URCK)=URCK-1) 
0807 
30 
CONTINUE 
0808 
C 
ADD ZONE JUST OUTSIDE K=KZ2 
0809 
ZHASCKZ2)=ZHASCKZ2)/2.0 
0810 
ZMASCKZ2+1)=ZMASCKZ2) 
0811 
V(KZ2+1)=V(KZ2) 
0812 
RHOZCKZ2+11=RHOZCKZ2) 
0813 
R(KZ2)=CCRCKZ2+1) •• J)-CJ.0.ZMASCKZ2+1). 
V(KZ2+1») •• CI.0/3.0) 
0814 
RZCKZ2)=CCRZ(KZ2+1) •• J)-CJ.0.ZHASeKZ2+1)/RHOZeKZ2+1») •• el.O/J.O) 
0815 


RRCKZ2)=RCKZ2)/RZCKZ2) 
0816 
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QRCKZ2+1)=RCKZ2+1)-RCKZ2) 
0817 
ORCKZ2)=RCKZ21-RCKZ2-1) 
0818 
R2CKl2)=RCKZ2) •• 2 
0819 
UCKZ2)=CC(UCKZ2-1)"2)+CUCKZ2+1) •• 21)/2.0).'0.S 
0820 
TRCKZ2+11=TR(KZ2) 
0821 
PCKl2+tl=P(KZ2) 
0822 
QCKZ2+11=Q(KZ2) 
0823 
TCKZ2+11=0.9tTCKl21 
0824 
T(KZ21=1.ltTCKZ21 
0825 
ECKZ2+11=O.8S.ECKZ21 
0826 
ECKZ2)=1.lS'ECKZ2) 
0827 
DIVFACKZ2+1)=DIVFACKZ2) 
0828 
URCKZ21=O.5tCURCKZ2-tl+URCKZ2+1» 
0829 


RETURN 
0830 


ENO 
0831 
• 
FORTRAN 
0832 
CFtCGSPO 
CGS PRINT OUT PACKAGE 
FIRE 1 
0833 
SU8ROUTINE CGSPO 
0834 
C 
STANDARD DIMENSION AND COMMON STATEMENTS 
AS IN MAIN PROGRAM 
0835 
45 
FORMAT(lH11139H AT COMPLETION OF MASTER CYCLE NUM8ER 
14.25X.15. 
0836 
45 
240X.3HSE=OPFIO.3116H TIME=IPEIO.3.21X.6HPOWER=lPEIO.3.2X.5HWATTS. 
0837 
45 
3t3X.tOHTIME STEP=lPEI0.3.20X.14HPASSES ENG Ea 141/ 
0838 
45 
43H 
K.2X.6HRADIUS.5X.8HPART VEL.2X.4HPRES.6X.3Ha/P.3X.3HTEV.4X.3HR 
0839 


45 
5HO.7X.3HE/G.7X.4HFOUT.5X.5HFI/FO.2X.4HDIVF.5X.4HTEMP.4X.5HGAMMA. 
0840 
45 
6X.4HPART.X.1HC.2X.3HTAU.5X.4HDTAUIICI4.1PE10.3.1P2EIO.2.0P2F6.2. 
0841 


4S 
71P3EIO.2.0PF6.2.1P2EIO.2.0PFS.2.0PF6.2.IAI.OPF7.2.1Al.lPE9.2)) 
0842 


46 
FORMATCtH IBH T PART=lPE12;4.9H FRACT 
.5X.6HDTHYD=IPE12.4.BH SEC 
0843 


46 
tS 
.5X.6HPOW03=lPEI2.4.BH WATTS 
.5X.5HTEFF=OPF6. O.BH DEG K 
.5X. 
0844 


46 
24HLZR=14/BH TYIELD=lPEI2.4·9H ERGS 
.5X.6HDTRAD=IPEI2.4.8H SEes 
0845 


46 
3 
.5X.6HPOW34=IPE12.4.BH W,\TTS 
.5X.5HTCOL=OPF6.0.BH DEG K 
.5X. 
0846 
46 
44HNCW=14/BH IN ENG=IPE12.4.9H ERGS 
.5X.6HDTMIN=IPEI2.4.8H SEes 
0847 


46 
5 
.5X.6HPOW45=IPEt2.4.BH WATTS 
.5X.5HWLMX=OPF6.0.BH ANG 
.5X. 
0848 


46 
64H~Dt=14/BH KN ENG=IPEI2.4.9H ERGS 
.5X.6HDTLST=IPEI2.4.BH SECS 
0849 


46 
7 
.5X.6HPOW57=IPEt2.4.BH WATTS 
.5X.5HTY4t=OPF7.2. 7HKT 
.5X. 
0850 


46 
B.4HE't1S=OPFS,4) 
0851 


36 
FORHATCBH TOT ES=IPEt2.4.9H ERGS 
5X.6HFBRAD=tPEI2.4.8H eM 
0852 


36 
25X.6HPOW71=lPEI2.4.BH WATTS 
.5X.5HTY71=OPF7.2. 7HKT 
.5X. 
0853 


36 
34HNSZ=14/BH E AMBT=IPE12.4.9H ERGS 
5X.6HSHRAD=IPEI2.4.BH CM 
0854 


36 
4 .5X.6HPOW47=IPE12.4.8H WATTS 
.5X.5HTY47=OPF7.2. 7HKT 
.5X. 
0855 


36 
S6HD15B9-121 
0856 


47 
FORMAT(lHlll. 
0857 


47 
43H 
K.2X.6HRADIUS.5X.BHPART VEL.2X.4HPRES.6X.3HQ/P.3X.3HTEV.4X.3HR 
0858 


47 
SHO.7X.3HE/G.7X.4HFOUT.5X.5HFI/FO.2X.4HDIVF.5X.4HTEMP.4X.5HGAHMA. 
0859 


47 
6X.4HPART.X.IHC.2X.3HTAU.SX.4HDTAUIICI4.1PE10.3.1P2E10.2.0P2F6.2. 
0860 


47 
71P3E10.2.0PF6.2.1P2E10.2.0PF5.2.0PFS.2.1Al.OPF7.2.1Al.lPE9.21) 
0861 
YIELD=EQUAL 
0862 
00 1000 K=4.LZ 
0863 
IFCTCK-l)-O.StTOOl 
1002.1002.1000 
0864 
1 002 
I(Z2= 0 
OB6~ 
CALL 0 I AGN$ 
... ~.-: 


1 auu 
CONT lNUEmf''l 
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C 


1059 


1060 
1061 
1062 
1064 
1063 
1065 
1065 
1067 
1068 
1164 
1165 
1165 
1167 
1168 
1264 
1265 
1265 
1267 
1268 
1364 
1365 
1365 
1367 
1368 
1464 
1465 
1465 
1467 
1468 
1564 
1565 
1565 


RSHOCK=O.O 
SAVE1=DT 
DT=CNCLZ) 
STX=PCLZ) 
POWER=TFLUX.l.0E-07 
TCOL=CTFLUX/( 
EMS 
.12.56.5.67E-05.(FBR •• 2))) •• 0.25 
TEFF=(TFLUX/(12.56.5.67E-OS·CFBR •• 2))) •• 0.25 
CALL PHOTOG 


ESTIMATE OF OBSERVED SPECTRAL DISTRIBUTION 
IF(FBR) 
1059.1059.1061 
P47=0.37.POWER 
047=047+P47.(TIME-TIMES)/4.186E+12 
041=0.0 
P03=0.0 
P34=0.0 
P45=0.0 
P57=0.0 
P71=0.0 
071=0.0 
GO TO 1600 
SAVE T = T( MCP) 
T( MCP) = TCOL 
IF( T(MCP)tCWLM 
-2.0E+08) 
1065.1065.1068 
CWLM=CWL( NCW. 1 ) 
FLM = EXPF((-3.41E-08tCWLM 
• T(MCP)+2.5E-25.((CWLM 
t T(MCP))tt3) 
l)tEXPF(-3.07E+13/((CWLM 
• T(MCP))t.I.B))) 
GO TO 1164 
FLM =7.6IE+22/((CWLM 
tT(M('P)) .. 3) 


IF( T(MCP)t3000.0-2.0E+08) 
1165.1165.1168 
F3 
= EXPF((-3.41E-08e3000.0t T(MCP)+2.5E-25t((3000.0t T(MCP))"3J 
1).EXPF(-3.07E+13/((3000.0t T(MCP))'.1.8))) 
GO TO 1264 
F3 
=7.6IE+22/((3000.0.T(MCP))t.3) 
IFC TCMCP)t4000.0-2.0E+08) 
1265.1265.1268 
F4 
= EXPF((-3.4IE-08t4000.0' T(MCP)+2.5E-25t((4000.0. T(MCP)) •• 3) 


1)'EXPF(-3.07E+13/((4000.0. T(MCP)).tl.8))) 
GO TO 1364 
F4 
=7.61E+22/(C4000.0.T(MCP)).t3) 
IF( TCMCP).5000.0-2.0E+08) 
1365.1365.1368 
F5 
= EXPF((-3.41E-08.5000.0. T(MCP)+2.5E-25t((5000.0. T(MCP)).t3) 
l)tEXPF(-3.07E+13/(C5000.0t T(MCP))t.I.8))) 
GO TO 1464 
F5 
=7.61E+22/((5000.0.T(MCP))tt3) 


IF( TCMCP)t7000.0-2.0E+08) 
1465.1465.1468 
F7 
= EXPF((-3.41E-08.7000.0. T(MCP)+2.5E-25.((7000.0t T(MCP))t.3) 
l)tEXPF(-3.07E+13/(C7000.0. TCMCP)) •• 1.8))) 
GO TO 1564 
F7 
=7.61E+22/((7000.0.TCMCP)) •• 31 
IF( T(MCP).10000.0-2.0E+08) 
1565.1565.1568 
FI= EXPF((-3.41E-08.10000.0tT(MCP)+2.5E-25.((10000.Ot T(MCP)).t3) 
1).EXPF(-3.07E+13/((10000.0. T(MCP)) •• 1.8))) 
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0868 
0869 
0870 
0871 
0872 
0873 
0874 
0875 
0876 
0877 
0878 
0879 
0880 
0881 
0882 
0883 
0884 
0885 
0886 
0887 
0888 
0889 
0890 
0891 
0892 
0893 
0894 
0895 
0896 
0897 
0898 
0899 
0900 
0901 
0902 
0903 
0904 
0905 
0906 
0907 
0908 
0909 
0910 
0911 
0912 
0913 
0914 
0915 
0916 
0917 
0918 


1567 
1568 
1580 


1581 
1582 
1583 
1584 
1585 
1586 
1587 
1588 
1589 
15S0 
15S1 


1600 


B 


488 


489 
490 


491 
492 


493 


GO TO 1580 
FI 
=7.6IE+22/CCI0000.0.TCMCP)) •• 3) 


CONT INUE 
TCMCP)=SAVET 
IFCCWLM-3000.0) 1583.1583.1582 
F3=FLM 
IFCCVLM-4000.0) 1585.1585.1584 
F4=FLM 
IFCCWLM-5000.0) 1587.1587.1586 
F5=FU1 
IFCCWLH-7000.0) 1589.1589.1588 
F7=FLM 
IFCCWLM-IOOOO.O) 1591.1591.1590 
FI=FlH 
. 
CONTINUE 
P03=POVER*CFLH-F3)/FLH 
P34=POVER*CF3-F4)/FlM 
P45=POVER*CF4-F5)/FLH 
P57=POVER*CF5-F7)/FLH 
P71=POWER.CF7-Fl)/FLH 
P47=POVER*CF4-F71/FLH 
071=071+P71*CTIME-TIMESl/4.186E+12 
047=047+P47*CTIME-TIHES)/4.186E+12 
041=047+071 


ENC 1 )=4. 189.CRC 1 )"31.CEC 1 )/VC 1 1 ) 
BNCI)=ENCI) 
CNC 1 )=3. 14159.CUC I ) .. 2hZHASC 1 1 
ONCI1=CNCI) 
OPPCl)=BNCI)+ONCI1 
00 490 
K=2.lZ 


ENCK1=4.18S*CCRCK) •• 3)-CRCK-l) •• 3)).CECK)/VCK)) 
CNCK)=3.1415S.CUCK) •• 2+UCK-I)·.2).ZHASCK) 
BNCK)=BNCK-l1+ENCK) 
ONCK1=ONCK-11+CNCK1 
THETACK)=212121212121 
IFCOCK)-STX) 
490.490.488 


STX=OCK) 
J3=K-3 
RSHOCK=RCK) 
OPPCK)=BNCK)+ONCK) 
EAMB=4.189.CRClZ) •• 3).EC99)/VC99) 
ESYSM=OPPClZ)-EAMB 
ETOTAL=ESYSM/4.1B6E+19 
PARTT=FLEX/YIElO 
IF(NMC) 
493.492.493 


YIElO=ESYSM 
EQUAL=YIELD 
CONT IHUE 
WRITE OUTPUT TAPE22.495.NR.NMC.TIME.EAHB.NCW.HCP.OPPCNCW).OPPCHCP) 
1.TRCHCV).TRCMCP).CK.ENCK).8NCK).CNCK).ONCK).OPPCK).ACK. 1).ZOCK. I). 
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0919 
0920 
0921 
0922 
0923 
0924 
0925 
0926 
0927 
0928 
0929 
0930 
0931 
0932 
0933 
0934 
0935 
0936 
0937 
0938 
0939 
0940 
0941 
0942 
0943 
0944 
0945 
0946 
0947 
0948 
0949 
0950 
0951 
0952 
0953 
0954 
0955 
0956 
0957 
0958 
0959 
0960 
0961 
0962 
0963 
0964 
0965 
0966 
0967 
0968 
0969 


495 
49S 
495 
49S 
-t95 
8494 
8 


828 
29 


496 
80 
81 
81 
81 
82 
92 
82 
83 
30 
30 
30 
85 
85 


-too 


2ZICK.I).TRCK).K=I.LZ) 
FORMATCIHI. 13H ENERGY CHECK.IOX.4HRUN=I-t.SX.4HNMC=14.SX.SHTIME=IPE 
210.3.SX.6HE AM8=IPEIO.31ISH NCW=14.SX.4HMCP=14.SX.9HDPPCNCW)=IPEI2 
3.4.SX. 9HDPPCMCP)=IPEI2.4.SX,8HTRCNCW)=IPEIO.3.SX,9HTRCMCP)=1PEtO 
4.311SH ZONE.3X.2HDE.IOX.4HSUME.6X.3HDKE.7X.SHSUMKE.SX.9HSUMTOTALE. 
SllX. IHA. II X. 2HZO. lOX. 2HZ I. lOX. 2HTRIIC IS. IPSEI2. 4, 3X. IP-tEI2. 4)) 


THETAC I 1=646464646464 
THETAC21=232323232323 
DO 29 
K=l.LZ 
HZ(K1=DTAUCK,11/l.S 
DMMCK1=OTAUCK.I1/CRCK1-RCK-l11 
FZCK1=4.28 E-09.(ECK1+P(K1.VCK1) 
OZCK)=TAUCK,I)/l.S 
OPCK)=QCKJ/PCK) 
FMCK)=TCK)/11606.5 
HPCK)=FOSCK).I.OE-07 
HPPCK)=A8SFCFISCK)/FOSCK)) 
GMCK)=1.0+CPCK).VCK1/ECK11 
URCK1=1.0/VCK1 
IFCPART(K)-4.0) 
29,28.29 


OTHCK)=3131313t3131 
CONTINUE 
< 
SHVEL=UCJ31.(1.0+GMCJ3))/2.0 
PSHK= C 2. O/C GMC J3) + 1.0) J .RH"ZC LZl. C SHVEL .. 2) 
ESHK=0.S-CC2.0-SHVEL/Cl.0+GMCJ3»))-·2) 
TRSHK=TCJ3J-ESHK/E(J3) 
ETOE=CS.67E-OS-CTRSHK.-4)/lRHOZCLZ).SHVEL»)+ECLZ) 
TTOE=ETOE-TCLZ)/ECLZ) 
POWSK=12.S6·CRCJ3+2)·-2).S.67E-OS-(TRSHK-.4J-FLM 
WRITE OUTPUT TAPE 6,496,NR.NMC,TIME,SHVEL.PSHK.TRSHK.TTOE 
PUNCH 496. NR.NMC.TIME,SHVEL.PSHK.TRSHK.TTOE 
FORMATC214,lPSEI2.4) 
IF(LZ-401 
30.30.81 
WRITE OUTPUT TAPE NTAPE.4S.NMC.LZ.ETOTAL.TIME.POWER.DT.NTI. 
I C K. RC K). UC K). PC K) • OPC K). FMC K). URC K). E C K). HP e K J • HPPC K) • 0 I VF AC K J • T C K 
2).GMCKJ.PARTCK).THETAeKJ.DZCKJ.DTH(K).HzeK).K=I.40) 
WRITE OUTPUT TAPE NTAPE.47, 
ICK,RCK),UCK),PCK),OPCK),FMCK),URCK),ECK),HPCK),HPPCK),OIVFACKJ.TCK 
2),GMCK),PARTCK),THETACK),OZ(K),DTH(K),HZCK),K=41,LZ) 
GO TO 85 
WRITE OUTPUT TAPE NTAPE,4S.NMC,LZ,ETOTAL,TIME.POWER.DT,NTI, 
ICK.RCK).ueK).PCK).DPCK).FM(K).URCK).ECK).HPCK).HPPCK).DIVFACK).T(K 
2).GMCK).PART(K).THETACK).DZCK).DTH(K),HZ(K).K=I.LZ) 
WRITE OUTPUT TAPE NTAPE.46.PARTT.RJ.P03,TEFF.LZR.FLEX,RX,P34.TCOL. 
2NCW.8NCLZ).DTMIN.P45.CWLCNCW.I).MCP.DNCLZ).SAVEI.PS7.Q41.EMS 
WRITE OUTPUT TAPE NTAPE.36.DPPCLZ).FBR.P71,Q71.NSZ,EAM8.RSHOCK.P47 
I.Q47.NR 
WCl.9)=RSHOCK 
DT=SAVEI 
IFCNMC-I) 
403.400.400 
IF C NMC -10) 
411 .411.403 
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0970 
0971 
0972 
0973 
0974 
0975 
0976 
0977 
0978 
0979 
0980 
0981 
0982 
0983 
0984 
0985 
0986 
0987 
0988 
0989 
0990 
0991 
0992 
0993 
0994 
0995 
0996 
0997 
0998 
0999 
1000 
1001 
1002 
1003 
1004 
1005 
1006 
1007 
1008 
1009 
1010 
1011 
1012 
1013 
1014 
1015 
1016 
lOt 7 
1018 
1019 
1020 


403 
CONTINUE 
1021 
C USER TAPE PRINT OUT 
1022 


~RITE OUTPUT TAPE 
32.300.N~.NHC.TIHE.LZ.P(99).TC99).RHOZC99).POWER 1023 
2.P47.FBR.CK.RCK).UCK).P(K).~RCK).T(K).FOSCK).FISCK).DHMCK1.O(K). 
1024 
3FZCK1.ECK).HPPPCK1.K=I.LZl 
1025 
300 
FORHAT(IHI.2I4. IPEI2.3.15.IP6EI2~3/CI4.IPI2EIO.3)) 
1026 
.WRITE OUTPUT TAPE 42.437.NR.NHC.TIHE.TCI).TCOL.TEFF.ETOTAL.PARTT. 
1027 
IBNCLZ).DNCLZ).ESYSH.EAHB.FLEX.041.Q47.071.P03.P34.P45.P57.P71.P47. 
1028 
2FBR.RSHOCK 
1029 


437 
FORMATCIH .2I4.IPIOEll.3/1PI2El1.4111) 
1030 
410 
IFCNMC) 
500.420.411 
1031 


411 
IFCTIHE-6.82E-OS) 
SOO.412.412 
1032 


412 
NSIX=NSIX+l 
1033 


413 
IF(NSIX-3) 
500.420.500 
1034 
420 
NSIX=O 
1035 


421 
CALL SIXDPO 
1036 
500 
RETURN 
1037 


END 
1038 
• 
FORTRAN 
1039 
CF1SIXD 
1040 


SUBROUTINE SIXDPO 
1041 
C 
STANDARD DIMENSION AND COMMON STATEMENTS 
AS IN MAIN PROGRAM 
1042 
FLAG=I.0E+31 
1043 
SOUNDS=SCRTFC1.4.P(99)/RHOZC99)) 
1044 
00 10 
K=I.LZ 
1045 
BN(K)=PCKJ/P(99) 
1046 
CNCK)=VC99)/VCK) 
1047 
DNCK)=TCKJ/T(99J 
1048 
ENCK1=UCKJ/SOUNDS 
1049 
C 
AVERAGE ZONE RADIUS IN FEET 
1050 
10 
DPCK)=CRCK)+RCK-I))/60.96 
1051 
WRITE TAPE 31.NR.TIME.LZ.CDP(K).BNCK).CNCK).DNCK1.ENCKJ.GMCK). 
1052 
IK=I.LZ).FLAG 
1053 
WRITE OUTPUT TAPE 41.50.NR.TIME.LZ.CDPCK).BNCK).CNCK).DNCK).ENCK). 
1054 
IGMCK).K=I.LZ) 
. 
lOSS 
SO 
FORMATC1HI/14.X.27H SIXDPLOT DIAGNOSTICS TIME=IPE10.3.5X.3HLZ=14 
1056 
SO 
211IC1P6EI6.4)) 
1057 
KZ6=KZ6+1 
1058 
20 
RETURN 
1059 


END 
1060 
• 
FORTRAN 
1061 
CFIPHOTOG 
PHOTOGRAPHIC BRIGHTNESS ROUTINE 
1062 


SUBROUTINE PHOTOG 
1063 
C 
STANDARD DIMENSION AND COMMON STATEMENTS 
AS IN MAIN PROGRAM 
1064 
IFCNMC) 
32.30.32 
1065 
30 
JTAPE=25 
1066 


31 
GO TO 40 
1067 


32 
IFCNHC- 3) 
33.33.30 
1068 


33 
JTAPE=6 
1069 
40 
CONTINUE 
1070 
C 
SOURCE FUNCTION EMCK) AND ABSORPTION COEFF BN(K) 
1071 


149 


DO 5 
K=I.LZ 
BNCK)=CVCK)/C 
16.S.(CI.293E-03.VCK)) •• I.I~)/(CCT(K)/11606.S) •• B.0)+2.6E-12.V(K») 
2+0. 180.(Cl.293E-03.VCK)) •• 1.91)/CCTCK)/11606.S) •• C3.S)) 
4+6.00E-OS.CCI.293E-03.V(K)) •• 2.0).CCTCK)/11606.S) •• 0.59) 
S+I.00E-09.CCI.293E-03.VCK)) •• 2.0).CCTCK)/11606.S) •• 4.00))) 
6.CI.0-EXPFC-3.2E+04/TCK»)) 
TEA=3.9S6E+04/TCK) 
IF( TEA-80. 0) 
I. 1.2 
2 
EI'1CK)=O.O 


3 
GO TO 5 


I 
EI'1CK)=I.OE+OS/CEXPF(3.9S6E+04/T(K))-1.0) 


5 
CONT INUE 


C OPTICAL DEPTHS 


DO 200 
K= I. LZ 
DI'1CK)=O.O 
DZeK)=O.O 
DPCK)=CRCK)+RCK-I))/2.0 
LS=K 
DO 100 
L=LS.LZ 
300 
IFCL-K) 
320.310.320 
310 
DZCL)=SQRTFCRCL) •• 2-DPCK).*2) 


311 
DI'1CL)=DZCL) 


31 3 
GO TO I 00 
320 
DZCL)=SQRTFCRCL)*.2-DPCK) •• 2)-DI'1CL-l) 


324 
DI'1(L)=DI'1CL-l)+DZCL) 
100 
HZCL)=EXPFC-BNCL).DZCL)/VCL)) 


C 
BRIGHTNESS INTEGRATION 
HPCLZ)=O.O 
LS=LZ-K+l 
DO 7 
J= 1. LS 


I'1=LZ-J 


7 
HP(I'1)=HPCI'1+1).HZCI'1+1)+Cl.0-HZCI'1+1)).EI'1CM+l) 
LS=K 
DO 6 
J=LS.LZ 


6 
HPCJJ=HPCJ-I).HZCJ)+CI.O-HZCJ)).EI'1CJ) 
HPPCK)=HP(LZ) 
(F(K-3) 
155.150.155 


150 
HNORI'1=HPPC3J 
IFCNMC) 
155.151.155 
151 
SCALE=HPP(3).10.0 


155 
CONTINUE 


199 
HPPPCK)=HPPCKJ/SCALE 
200 
HPPCK)=HPP(KJ/HNORI'1 
210 
WRITE OUTPUT TAPE 
JTAPE.201.NMC.TII'1E.CK.DPCK).HPPPCK).HPP(K).K=I. 
210 
ILZ) 
201 
FORI'1ATCIHI.6HPHOTOG. (4.IOX.5HTIME=IPEIO.31ISH ZONE.3X.6HI'1EAN R.4X. 


201 
25HABS B.SX.5HREL BIICI6.1P3EIO.3)) 
280 
RETURN 
END 
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1072 
1073 
1074 
1075 
1076 
1077 
1078 
1079 
1080 
1081 
1082 
1083 
1084 
1085 
1086 
1087 
1088 
1089 
1090 
1091 
1092 
1093 
1094 
1095 
1096 
1097 
1098 
1099 
I 100 
1 101 
1102 
1103 
1104 
1105 
1106 
1107 
1 t 08 
1109 
1110 
1 I I I 
1 1 12 
11 13 
1114 
1115 
1116 
1117 
111 B 
II 19 
1120 
1121 
1122 


• 
FORTRAN 
1123 


CFIDIAG 
1124 


SUBROUTINE OIAGNS 
1125 


C 
STANDARD DIMENSION AND COMMON STATEMENTS 
AS IN MAIN PROGRAM 
1126 
35 
FORMATCIHI.4HNMC=I4.5X.4HNTI=I4.5X.4HKZ2=15.5X.5HZONE=14.5X.4HRUN= 
1127 


214/) 
1128 


M=KZ5 
1129 
IF(KZ2) 
5.10.15 
1130 


5 
WRITE OUTPUT TAPE S.S.KZ2 
1131 


S 
FORMATC1HI.2IHOIAG CALLED FROM MAINI5) 
1132 


GO TO 20 
1133 


10 
WRITE OUTPUT TAPE S.ll.KZ2 
1134 


II 
FORMATCIHI.22HDIAG CALLED FROM CGSPOI5) 
1135 


GO TO 20 
1136 
15 
WRITE OUTPUT TAPE S.IS.KZ2 
1137 
16 
FORMAT(IHI.22HOIAG CALLEO FROM HYOROI5) 
1138 


20 
WRITE OUTPUT TAPE S.21.NMC.NTI.LZ.LZR.MCP.MCL.TIME.DT.CS.CR.BR.RS. 
1139 


20 
l(K.VCKl.VRCKl.ECK).ERCKl.RESDUECKl.PACKl.DIVFACKl.DIVFR(Kl.K=I.LZ) 
1140 


21 
FORMATCIH .4HNMC=I4.2X.4HNTI=14.2X.3HLZ=14.2X.4HLZR=14.2X.4HMCP=14 
1141 
21 
1.2X.4HMCL=14/SH TIME=IPEI2.4.2X.3HDT=IPEI2.4.2X.3HCS=OPF5.2.2X. 
1142 
21 
23HCR=OPF5.2.2X.3H8R=OPF5.2.2X.3HRS=OPF5.2113H 
K.6X.IHV.13X.2HVR. 
1143 
21 
312X.IHE.13X.2HER.12X.SHRESDUE.8X.2HPA.12X.5HOIVFA.9X.5HOIVFRII 
1144 
21 
4 ( I 4. 1 P8E 1 4 . 7 ) ) 
1 I 45 
22 
WRITE OUTPUT TAPE 6.23.NMC.NTI.OTMIN.YIELO.SCALE.TFLUX.FBR.TIMEW. 
1146 
22 
I(K.TCK).OTMCKl.OETCK).OPTCK).8NCKl.CNCK).ONCK).ENCK).K=I.LZ) 
1147 
23 
FORMATCIHl.4HNMC=14.2x.4HNTI=I4.2X.SHOTMIN=IPE12.4.2X.6HYIELO=IPEI 
1148 
23 
12.4.2X.SHSCALE=lPE12.4.2X/7H TFLUX=lPEI2.4.2X.4HFBR=IPEI2.4.2X.6HT 
1149 
23 
2IMEW=IPEI2.4113H 
K.SX.1HT.13X.3HOTM.IIX.3HOET.IIX.3HOPT.I1X.2HBN. 
1150 


23 
312X.2HCN.12X.2HON.12X.2HENIICI4.1PBEI4.7)) 
1151 


WRITE OUTPUT TAPE 6.35.NMC.NTI.KZ2.M.NR 
1152 
24 
WRITE OUTPUT TAPE 6.25.CK.FOS(K).FISCKl.BBCK.l).BCCK.I).RCK). 
1153 
24 
lZ0CK.1).A(K.ll.ZI(K.l ).K=I.LZ) 
1154 


25 
FORMAT(IH 13H 
K.6X.2HFO.12X.2HFI.12X.2HBB.12X.2HBC.12X.IHR.13X. 
1155 
25 
12HZO.12X.IHA.13X.2HZI/CI4.1P8E14.7)) 
1156 


WRITE OUTPUT TAPE 6.35.NMC.NTI.KZ2.M.NR 
1157 
26 
WRITE OUTPUT TAPE S.27.CK.SCK.l).WCK.l).OTAUCK.ll.R2(K). 
ZMASCKJ. 
1158 
26 
ITAUCK.I).CWLCK.ll.TRCKl.K=I.LZl 
1159 
27 
FORMATCIH 13H 
K.SX.lHS.13X.lHW. 13X.4HOTAU. 10X.2HR2.12X.4HZMAS. lOX 
1160 


27 
1.3HTAU.l1X.3HCWL.IIX.2HTR/CI4.1P8EI4.7)) 
1161 
WRITE OUTPUT TAPE S.35.NMC.NTI.KZ2.M.NR 
1162 
28 
WRITE OUTPUT TAPE S.29.CK.OMMCK).OMCK).OZ(K).OPCK).OPPCK).QCK). 
1163 
28 
lFOCK.2).FICK.2l.K=I.LZ) 
1164 
29 
FORMATCIH 13H 
K.SX.3HOMM.IIX.2HOM.12X.2HOZ.12X.2HOP.12X.3HOPP.IIX 
1165 
29 
1. IHQ.13X.5HFOTJZ.9X.5HFITJZ/CI4.IP8EI4.7)l 
1166 
WRITE OUTPUT TAPE S.35.NMC.NTI.KZ2.M.NR 
1167 


30 
WRITE OUTPUT TAPE 6.31.(K.FMM(K).FMCK).FZCK).HZ(K).FPCK).HP(K). 
1168 


30 
IHPP(K).HPPPCK).K=I.LZ) 
1169 
31 
FORMAT(IH 13H 
K.6X.3HFMM. IIX.2HFM. 12X.2HFZ. 12X.2HHZ. 12x.2HFP. 12X. 
1170 


31 
12HHP. 12X.3HHPP.llX.4HHPPP/CI4. IP8E14.7)) 
1171 
WRITE OUTPUT TAPE 6.35.NMC.NTI.KZ2.M.NR 
1172 
32 
WRITE OUTPUT TAPE 6.33.TDN2.TD02.TIMES.Q47.Q71.FLIX.FLOX.FLEX. 
1173 


lSI 


32 
lRH.RX.RJ. CK.FOCK.3).DTRClO.RZClO.PClO.DECK).RRClO.RHOZ(lO. UCK). 
1174 
32 
2K=1.100) 
1175 
33 
FORHATCIH .2X.SHTDN2=IPEI2.4.2X.SHTD02=IPEI2.4.2X.6HTIMES=IPEI2.4. 
1176 
33 
12X.4HQ47=IPEI2.4.2X.4HQ71=lPE12.~.2X/6H FL1X=IPE12.4.2X.5HFLOX=IPE 
1177 
33 
212.4.2X.SHFLEX=IPE12.4.2X.3HRH=IPE12.4.2X.3HRX=IPEI2.4.2X.3HRJ=IPE 
1178 
33 
312.4113H 
K.6X.3HDTH.1IX.3HDTR.11X.2HRZ.12X.1HP.13X.2HDE.12X.2HRR. 
1179 
33 
412X.4HRHOZ.I0X.2HU ICI4.1P8E14.7» 
1180 


RETURN 
1181 
END 
1182 


- 
FORTRAN 
1183 


CFIFLUX 
1184 


SUBROUTINE FLUXS 
1185 


C 
STANDARD DIMENSION AND COHMON STATEMENTS 
AS IN MAIN PROGRAM 
1186 


500 
LZMI=LZ-l 
1187 


00 530 
N=I.NSZ 
1188 


C 
BOUNDARY CONDITION 
NO INWARD FLUX AT OUTER BOUNDARY 
1189 


FICLZ.N)=O.O 
1190 
510 
DO 512 
J=I.LZMl 
1191 
511 
K=LZ-J 
1192 
512 
FICK.N)=FICK+l.N)-ZICK+I.N)+BBCK.I)-CACK+I.N»-2.0.SCK.N)-WCK+ 
1193 
512 
I I • N) 
1 194 


520 
FOCI.N)=FICI.N).Z(CI.N)+BBC1.2).CACI.N»+2.0.SCI.N)*WCI.N) 
1195 
522 
00 525 
K=2.LZ 
1196 
525 
FOCK.N)=FOCK-I.N).CR2CK-I)/R2CK».ZOCK.N) 
1197 
525 
2+BBCK.2).CACK.N»+2.0.S(K.N).WCK.N) 
1198 
525 
3+CI.0-CR2CK-l)/R2CK»).FICK.N).ZICK.N) 
1199 


530 
CONTINUE 
1200 


549 
00555 
K=I.LZ 
1201 
FISCK)=O.O 
1202 
FOSCKJ=O.O 
1203 
00 550 
N=I.NSZ 
1204 
FISCK)=FISCK)+FICK.N) 
1205 
550 
FOSCK) =FOSCK) +FOCK. N) 
1206 


551 
tF(K-l) 
552.552.554 
1207 


552 
01=0.0 
1208 


553 
GO TO 555 
1209 


554 
01=1.0 
1210 


555 
0IVFA(K)=CR2CK)/ZMAS(K»*CFOSCK)-FISCK» 
1211 
555 
I-Ol-(R2CK-l)/ZHASCK».CFOSCK-l)-FISCK-l» 
1212 
541 
FLOX=12.56-R2CLZ-3)-FOSCLZ-3) 
1213 


590 
RETURN 
1214 
591 
END 
1215 
1216 


• 
FORTRAN 
1217 


CFJSTE 
1218 


SUBROUTINE STATE 
1219 


C 
STANDARD DIMENSION AND COMMON STATEMENTS 
AS IN MAIN PROGRAM 
1220 
1221 


C 
STEFAN =2.S14E-09 TO INCLUDE RADIATION ENERGY AND PRESSURE 
1222 


C STEFAN =0.00 TO DELETE RADIATION ENERGY AND PRESSURE 
1223 


STEFAN=O.O 
1224 


152 


100 
DO 126 
K=I.LZ 
1225 


101 
ETACK)=1.0/Cl.293E-03' VCK» 
1226 


102 
TseKl=TseK)'I.OE-04 
1227 


103 
XSCK1=1.0E+04.TSCK)/( ETACK)'.0.086) 
1228 


104 
GNUeK)=LOGFeXSeK)/2000. O)/LOGFCI250. 0) 
1229 


lOS 
IF(XS(K)-2000.0) 106.106.109 
1230 
106 
PART(K)=1.0 
1231 


107 
GO TO 1 12 
1232 


109 
IFCXS(K)-2.5E+06) 
109.111.111 
1233 


109 
PART(K)=1.0+15.4'( GNUCK)"3)'C4.0-3.0' GNU(K» 
1234 


110 
GO TO 112 
1235 


III 
PART(K)=16.4 
1236 
112 
PCK)= (2.881.(TS(K)1 VCK»' PART(K)+STEFAN.(TS(K)"4».1.0E+10 
1237 
113 
Y=0.0794/C2.881.TSCK)' PARTCK» 
1238 
114 
UZ=I.0+C27.0'Y+3.0)/(5.0'Y+I.0)+961.0'Cl.0-Y)'V/(3000.0,(v'.2)+1.0 
1239 
114 
1)+C2356.0.Cl.0-Y)'V)/C9.0E+04,CY"2)+1.0)+C 41000.0'(1.0-V).V)1 
1240 


114 
2CI2.0E+06'CV"2)+1.0)+C2.15E+05'Cl.0-V).Y)/C1.5E+18'(V"4)+1.0) 
1241 
115 
UT=C24.0,CY"2)+4.0E-I0)/C4 .• (V •• 2)+I.E-I0)-CO.0970.CV"2)'C1.0-V) 
1242 
115 
1)/C2.0E-06+Vu3)+C4.19E-05HVu3).CI.0-V»/(I.14E-11+vu6) 
1243 
116 
UC=UZ-0.09'(UZ-UT)'LOGFC ETACK» 
1244 
117 
E(K)= (0.03920.(UC-l.0)/V+3.0'STEFAN.(TS(K) •• 4)' V(K».1.0E+10 
1245 
126 
TS(K)=TS(K).1.0E+04 
1246 
TS(2)=TS(2)'1.0E-04 
1247 


TS(I)=TSCI).1.0E-04 
1248 


P(I)=CO.538E+10.TS(I)"1.5)/Y(I) 
1249 


ECI)=2.53E+l0'(((Y(I)"0.25)+10.0)/(3.0'Y(1)"0.5+10.0». 
1250 
I(CV(I)"0.25)+0.13).(TS(I)"1.5+.02722/V(I)"1.5) 
1251 


P(2)=36.18.TS(2).1.0E+I0.(920.0+TS(2).'2)/(V(2).(TS(2)"2+1.08E+04 
1252 


I» 
1253 
E(2)=((649.0+TS(2)"2)/(100.0+TS(2») 
1254 


I.CC82.7.TS(2)'I.OE+I0)/(116.0/CV(2)'.0.25)+TS(2)'(1.0+0.12/(V(2) •• 
1255 


20.25»» 
1256 


TS(I)=TS(I).1.0E+04 
1257 


TS(2)=TSC2)'1.0E+04 
1258 
176 
RETURN 
1259 
177 
END 
1260 
• 
FORTRAN 
1261 


CFIHYDO 
1262 


SUBROUTINE HVDRe 
1263 
C 
STANDARD DIMENSION AND COMMON STATEMENTS 
AS IN MAIN PROGRAM 
1264 
1265 


KHR=O 
1266 
KRS=1 
1267 


KZ4=0 
1268 
C 
TEST FOR IMPROPER CONVERGENCE 
1269 


DO 195 
K=I.LZ 
1270 
IFCDTMCK» 
192.195.191 
1271 
191 
IFCTCK)-2.0E+05) 
188.188.189 
1272 
198 
IFCABSFCDTMCK»-4.99.TCK» 
195.195.193 
1273 
199 
IFCABSFCDTM(K»-0.80.TCK» 
195.195.193 
1274 
192 
IFCABSFCDTMCK»-O.25.TCK» 
195.195.193 
1275 
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193 
KZ2=193 
1276 


RQX=RQX.0.2 
1277 


RH=RQX.I.3 
1278 


CALL OIAGNS 
1279 
194 
GO TO 198 
1280 


195 
CONTINUE 
1281 


198 
IF(RQX) 
199.199.200 
1282 
199 
RQX=DTMIN 
1283 


C RICHTMYER 
VON NEUMANN HYDRODYNAMIC SCHEME 
1284 
C 
ADVANCE VELOCITY - - CHOOSE TIME STEP 
1285 


200 
DO 207 
K=I.LZ 
1286 


201 
"CK.31=UCKl 
1287 
202 
URCK1=UCKl 
J288 


203 
UCK)= UCKl-2.0.CCR2CKl 
.DT1/(ZMASCK1+ZMASCK+ll»).CPCK+1)+Q(K+1J 
1289 


2041-PCKJ-QCK» 
1290 
206 
IFCUCK)-1.0E+11) 207.210.210 
1291 


207 
CONTINUE 
1292 


208 
UCLZ1=0.0 
1293 


209 
GO TO 217 
1294 
210 
KZ2=210 
1295 


CALL D I AGNS 
1296 
213 
GO TO 52 
1297 


C HYDRODYNAMIC TIME STEP - - COURANT CRITERION 
1298 
217 
RM=RQX-1.3 
1299 
218 
IFCRM-0.070.TIMEJ 
222.222.219 
1300 


219 
RM=0.070-TIME 
1301 


222 
DO 227· K=I.LZ 
1302 


224 
DTHCKJ=DRCK)/CRS.SQRTFCP(KJ/(VCK).CRHOZCK) •• 2)))) 
1303 
DTH(K)=DTH(K)-C(RZCK1/RCK»).-2) 
1304 
FOCK.3)=DTHCK) 
1305 


225 
IHRM-DTHCK1) 227.226.226 
1306 


226 
RM=DTHCKl 
1307 
KZ4=K 
1308 


227 
.CONTINUE 
1309 


228 
RJ=RH 
1310 
1311 


C 
RADIATION LIMITED TIME STEP 
1312 
229 
RX=RM 
1313 


230 
DO 242 
K=KRS.LZ 
1314 
IHK-l) 
233.233.223 
1315 
223 
IFCRCK1-FBR) 
231.231.233 
1316 
231 
AR=BR 
1317 
232 
GO TO 234 
1318 
233 
AR=BR.2.0 
1319 
234 
CONTINUE 
1320 


235 
DTRCK)=ABSFCAR-ECK)/DIVFACK») 
1321 
238 
IFCDTRCK1-OTMIN 
1 242.242.239 
1322 
239 
IHOTRCK1-RXl 
240.240.242 
1323 


240 
RX=DTRCK) 
1324 
241 
KZ4=K 
1325 


242 
CONTINUE 
1326 
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243 
RH=RX 
244 
IFCRH-2.0eDT) 
246.246.245 
245 
Rt1=2.0eDT 


246 
IFCRt1-DTHIN) 
247.248.248 
247 
Rt1=DTHIN 
248 
RQX=RH 


249 
IFC(TIHE+RM1-TIHEW1 
252.252.250 
250 
IFCNMC- 41 
252.252.251 
251 
RH=TIMEW-TIME 


252 
DO 255 
L=I.LZMI 


253 
UACL1=CUCLl/2.01+CURCLl/2.0)+CUCL1.RMl/C2.0-0Tl-CURCL).RMl/C2.-0T) 


254 
URCL1=CUCLl/2.0)+CURCLl/2.0l-CUCL).RM)/C2.0.DT)+CURCL).RM)/(2.-0T) 


255 
UC U =UA( U 


256 
DT=RM 
C 
ADVANCE RADIUS ALL L 


260 
DO 263 
L=I.LZ 
WCL.5)=RCU 


261 
RCLl=RCL)+OT.UCL) 


262 
RRCL)=RCL)/RZCL) 


263 
R2CL)=RCL) •• 2 
00 269 
K=2.LZ 
IFCCRCK)-RCK-l».,.0.15) 
264.264.269 
264 
KHR=KHR+I 
IFCKHR-3) 
265.265.270 


265 
KZ2=266 


266 
CALL OIAGNS 
00 268 
L=I.LZ 
UCU=WCL.3) 
RCU=WCL.5) 


268 
CONTINUE 
RQX=DhO.2 
GO TO 200 


269 
CONTINUE 
C 
ADVANCE SPECIFIC VOLUME 
270 
V(1)=CCRCI)--3)/CRZ(I) •• 3»/RHOZC1) 


271 
00 
277 
K=2.LZ 


273 
VCK)=(I.O/RHOZCK))-CCCRCK)/ORCI)) •• 3)-CCRCK-I)/ORCIll- -3)1 


274 
ICCCRZCK)/ORCI))-.3l-CCRZCK-I)/ORCI)) •• 3ll 
275 
IFCVCK)-I.OE+20l 
276.279.279 
276 
IFCVCK)) 
279.279.277 
277 
CONTINUE 


278 
GO TO 283 
279 
KZ2=279 
CALL DIAGNS 


282 
GO TO 52 


C 
ADVANCE ARTIFICAL VISCOSITY 


283 
IFCNMC-NQS1 
284.284.286 
284 
VD=CS 


285 
GO TO 287 
286 
VD=CR 


287 
00 295 
K=I.LZ 
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1327 
1328 
1329 
1330 
1331 
1332 
1333 


, 1334 
1335 
1336 
1337 
1338 
1339 
1340 
1341 
1342 
1343 
1344 
1345 
1346 
1347 
1348 
1349 
1350 
1351 
1352 
1353 
1354 
1355 
1356 
1357 
1358 
1359 
1360 
1361 
1362 
1363 
1364 
1365 
1366 
1367 
1368 
1369 
1370 
1371 
1372 
1373 
1374 
1375 
1376 
1377 


IFeVCK1-VRCKll 293.2942.2942 
293 
Q(Kl=CVD.CCZMASCK1/RCK1 •• 2) •• 21/(VCK)+VRCK1».CCCVCK1-VRCK»/DT) •• 
12) 
294 
IFCQ(K)-1.0E+22) 2941.298.298 
2941 
IFC 
QCK» 
2942.295.295 
2942 
Q(K)=O.O 


295 
CON T I NUE 
296 
QCLZ)=O.O 


297 
RETURN 
298 
KZ2=298 
CALL DIAGNS 


52 
WRITE OUTPUT TAPE 6.43.Tel) 


43 
FORMAT(IH .15HEXIT FROM HYDRO.5X.5HTCI1=IPEI3.6) 
T( I ) = 0.000000 
1500 
GO TO 297 
END 


• 
FORTRAN 
CF1COEF 
SUBROUTINE 
COEFF 
C 
STANDARD DIMENSION AND COMMON STATEMENTS 
AS IN MAIN PROGRAM 
700 
LZR=LZR 
C SET FLUX DERIVATIVES NEAR OUTER BOUNDARY OF RADIATIVE REGION 
701 
FP(LZR)=O.O 
702 
HPCLZR)=O.O 
703 
HPPCLZR)=O.O 
704 
HPPPCLZR)=O.O 
705 
HPP(LZR-Il=O.O 
706 
HPPP(LZR-I)=O.O 
707 
HPPP(LZR-21=0.0 
DO 769 
J=I.LZR 
C 
SET RIPPLE ZONE PARAMETERS 
DO 
514 
N=I.NSZ 
SAVEI=BBeJ.N) 
SAVE2=BB( J -I. N) 
SAVE3=SCJ.N) 
SAVE4=SeJ-I.Nl 
SAVE5=DTAUCJ.N) 
SAVE6=AeJ.N) 
SAVE7=WCJ.N) 
SAVE8=Z IC J. N) 
SAVE9=ZI(J-I.Nl 
SAVE10=ZOCJ.Nl 
SAVEII=ZOCJ+I.N) 
SAVEI2=BCCJ.N) 
SA VE I 4 = T( J ) 
SAVEI3=TSeJl 
SAVE 1 5=CWL( J. 1 ) 
SAVEI6=ZOeJ+2.N) 
SAVEI7=BBCJ.2) 
SAVEIB=BBCJ-l.2l 
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1378 
1379 
l3BO 
1381 
1382 
1383 
1384 
1385 
:386 
1387 
1388 
1389 
1390 
1391 
1392 
1393 
1394 
1395 
1396 
1397 
1398 
1399 
1400 
1401 
1402 
1403 
1404 
1405 
1406 
1407 
1408 
1409 
1410 
I 4\ I I 
1412 
1413 
1414 
1415 
1416 
1417 
1418 
1419 
1420 
1421 
1422 
1423 
1424 
1425 
1426 
1427 
1428 


311 
BBCJ.Nl=CCCCCRCJ+I)-RCJl).TPCJl •• 4)+CCReJl-RCJ-I».T eJ+I) •• 4»1 
1429 
311 
lCRCJ+I)-ReJ-I»».S.67E-OS 
1430 
BBCJ-l.Nl=CCCCCRCJ)-RCJ-I».TCJ-I) •• 4)+CCRCJ-I)-RCJ-2)).TPeJ) •• 4» 
1431 


l/CReJ)-RCJ-2»l).S.67E-OS 
1432 


B8CJ-l.2)=BBCJ-l.1) 
1<433 
BBCJ.2)=BBeJ.I) 
1434 
312 
BCCJ.N)=S.67E-OS.eTPCJ) •• 4) 
1435 


TSCJ)=TPCJ)/11606.5 
1436 


TeJ)=TPCJ) 
1437 


ETACJ)=1.293E-03.VCJ) 
1438 


ALPHA=9.0+EXPFC-TIME/I.OE-03.CETACJ) •• 2» 
1439 


8ETA=CC6.0E-18.VCJ» •• 0.5)'CTSCJ) •• 0.1) 
1440 
300 
ZICJ.I)=EXPFC-I.5.CRCJ)-RCJ-I» 
ICCETACJ) •• 1.5)/CCTSCJ) •• ALPHA) 
1441 


300 
2+BETA) 
+0.2.CETACJ) •• 1.91)/CTSCJ) •• 2.73) 
1442 


300 
3+0.023.CETACJ) •• 1.8).CTSCJ) •• 0.25)+I.OE-07.CETACJ) •• 2.0)1 
1443 


300 
4CCTSCJ) •• C-5.0»+4.0E-14.VCJ»1) 
1444 


ZICJ.I)=ZICJ.l1.0.99999998 
1445 


Z I C J • 2) = Z I C J. 1 1 
1446 
HMFP=CETACJ1 •• l.51/CCTSCJ) •• ALPHA)+BETAl 
1447 


1+0.2.CETACJ) •• 1.91)/CTSCJ) •• 2.73) 
1448 


2+0.023.CETACJ) •• 1.8).CTSCJ) •• 0.25) 
1449 


3+I.OE-07.CETACJ) •• 2.0).CTSCJ) •• S.00) 
1450 
HMFF=9.3E-06.CETACJ) •• 2.0).CTSeJ) •• 3.5) 
1451 


IFCTSCJ)-10.0) 
301.298.298 
1452 
298 
IFCHMFF-HMFP) 
299.301.301 
1453 
299 
HMFP=HMFF 
1454 


301 
OTAUeJ.I)=1.5.CRCJ)-ReJ-I»/HMFP 
1455 


302 
ACJ.N)=EXPFe-OTAUeJ.N» 
1456 


305 
TSeJ)=TPCJ) 
1457 


506 
SCJ.N)=CBCCJ.N)-BCCJ+I.N»/COTAU(J.N)+OTAUCJ+I.N» 
1458 


SeJ-I.N)=CBCCJ-I.N)-BCeJ.N»/eOTAUCJ-I.N)+OTAUeJ.N» 
1459 


507 
WCJ.N)=I.O-ACJ.N)-ACJ.N).OTAUeJ.N) 
1460 
508 
IFeWeJ.N)-I.OE-04) 
509.514.514 
1461 


509 
WCJ.N)=0.5.COTAUCJ.N) •• 2) 
1462 
514 
CONTINUE 
1463 
1464 


C 
SET CUTOFF WAVELENGTHS 
1465 


540 
I =J 
1466 
541 
IFCTCI1-TON2) 
546.542.542 
1467 
1468 


C 
OPACITY DUE TO ATOMIC SPECIES 
1469 
542 
CWLel.l)= 700.0.EXPFC-0.36.CTCI)-11606.5)/11606.5) 
1470 
IFCCWLCI.I)-275.0) 
545.560.560 
1471 
545 
CWLCI.Il=275.0 
1472 


GO TO 560 
1473 


5<46 
IFCTCI1-T002) 
551.547.547 
1474 
1475 


C OPACITY DUE TO N2 MOLECULE 
1476 
547 
TCN2=5.0E-07.1.293E-03.VCI) 
1477 


CWL( I • I ) = 1140. O. C C C R C I ) -RC 1-1 )) I C V C I ).1 • 293E - 03))" O. II ) 
1478 


1.CI.O-EXPFC-TIME/TCN2» 
1479 
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IFCC~LCI.I)-1000.0) 550.560.560 
1480 
550 
CVlCI.l)=1000.0 
1481 


GO TO 560 
1482 
1483 


C 
OPACITY DUE TO 02 MOLECULE 
1484 
SSt 
TC=42S0.0-271.0.LOGFCI.293E-03. YCI)) 
1485 


TC02=3.OE-07.1.293E-03.YCI) 
1486 
553 
IFCT(I)-TC)SS4.SS4.SS6 
1487 
554 
C~L(I.I)=IS00.0+TS(I).CO.163+0.074S.LOGF(DR(I).RHOZCI)11.293E-03» 
1488 
554 
t.Cl.0-EXPFC-TIHE/TC02) 
1489 
555 
GO TO 557 
1490 


556 
TSB=TSCI).C.0647.LOGFCRCI)-RCI-I»-0.2S-.109.LOGFC1.293E-03.Y(I») 
1491 


C~LCI. 1)=3S00.0+TSB-C2000.0+TSB).EXPFC-TIME/TC02) 
1492 
557 
IFCCWLCI.Il-ISOO.O) 558.560.560 
1493 
558 
CWLCI.I)=IS00.0 
1494 


560 
CONTINUE 
1495 
1496 


C 
CALCULATE ZO VALUES 
1497 


ZO C 1 • 1 ) = Z I C 1 • 1 1 
1 498 


604 
JSTR=J+2 
1499 


KD=O 
1500 
60S 
DO 670 
L=J. JSTR 
1501 
TRCL1=TCLl 
1502 
IF(TCLl-8.0E+04) 
607.668.668 
1503 


607 
KD=KD+I 
1504 


IFCKD-IJ 
608.608.609 
1505 


608 
HS=L-l 
1506 


609 
IFCOTAU(L.IJ-I.OJ 
610.668.668 
1507 


C 
SELECT TRCL)= TEHPERATURE OF RADIATION TRAVERSING ZONE L 
1508 
610 
LIH=L-l 
1509 
61 I 
T AUSUH = o. 0 
I 5 1 0 


6 1 5 
00 630 NN = I • LI M 
I 51 1 
616 
H=L-NN 
1512 
617 
TAUSUM=TAUSUM+DTAUCH.I) 
1513 


620 
IFCTAUSUM-0.7J 
630.627.627 
1514 
627 
TRCL)=TCMJ 
1515 
628 
IFCH-2) 
630.630.638 
1516 


630 
CONTINUE 
1517 


IFCTCLJ-2500.0) 
631.632.632 
1518 
631 
TRCL)=TRCL-I) 
1519 


GO TO 637 
1520 
632 
EPTSUH=O.O 
1521 


EPSUH=O.O 
1522 


00 635 
1=3.LlH 
1523 


EPTSUH=EPTSUM+TCI).DTAUCI.IJ 
1524 


EPSUH=EPSUH+DTAUCI.l) 
1525 


IFCT(I)-2500.0) 
636.636.635 
1526 
635 
CONTINUE 
1527 
636 
TRCL)=EPTSUH/EPSUH 
1528 
637 
H=2 
1529 
638 
I Z=H+ 1 
1530 
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639 
DO 646 
I=IZ.L 
1531 


640 
IF(TR(L).CWLCI.I)-2.0E+08) 
641.641.644 
1532 
641 
ZI(I.2)= EXPF((-3.41E-08.CWLCI.l).TRC L )+2.5E-25.C(CWLCI.l).TRCL 
1533 
641 
2» •• 3».EXPFC-3.07E+13/CCCWLCI.t).TR( L » •• 1.8») 
1534 
642 
GO TO 646 
1535 


644 
ZICI.2J=7.61E+22/CCCWLel.l).TRCL » •• 3) 
1536 
646 
CONTINUE 
~537 
ZMAX=ZIeM+l.2) 
1538 


ZOeLo1)=ZI(M+I.2) 
1539 
647 
IFCL-M-l) 668.663.650 
:540 
C 
CORREC TI ON FOR I MTERMED I ATE ZUNES. 
541 


650 
IZ=M+2 
1542 
651 
CWLTST=CWLCM+I.l) 
1543 
652 
DO 660 
NZ=IZ.L 
1544 
653 
IFCCWLCNZ.l)-CWLTST) 
654.654.657 
1545 
654 
IFCNZ-L) 
660.655.668 
1546 
655 
ZOCL.l)=1.0 
.1547 
656 
GO TO 660 
1548 
657 
ZOCL,I)=Z)CNZ.2)/ZMAX 
1549 
658 
CWLTST=CWLCNZ.l) 
1550 


659 
ZMAX=ZICNZ.2) 
1551 


660 
CONTINUE 
1552 
663 
)FCZOCL. I)-ZICL. I» 
670.668.668 
1553 
668 
ZOCL.I)=ZI(L.l) 
1554 


670 
CONTINUE 
1555 
671 
IFCDTAUeJ.l)-1.0E-04) 
672.674.674 
1556 
672 
A(J.l )=DTAUeJ.l) 
1557 
673 
GO TO 675 
1558 
674 
ACJol )=1. O-AeJ.l) 
1559 
675 
IFeSCJ.l)-1.69E+38) 
8675)8675.677 
1560 


8675 
1F((BBCJ.l).A(J+l.l)-2.0.S~J.l).weJ+I.I))-BCCJ+I.I).ACJ+1.1).0.01) 1561 
8675 1 677.680.680 
1562 
677 
S(J.l)=O.O 
1563 
678 
BBCJ.I)=BCeJ+l. I) 
1564 
679 
BBeJ.2)=BCCJ.I) 
1565 


680 
IFeS(J-I.I)-1.69E+38) 
8680.8680.681 
1566 


8680 
I F C ( BB e J - 1 • I ). A e J. 1 ) - 2. O. S e J - 1 • 1 ). we J. 1 )) -BC e J. 1 ) It A e J • 1 ) • O. 0 I ) 
1 567 


8680 I 
681.684.684 
1568 
681 
S (J -1 • 1 ) = O. 0 
1569 
682 
BBeJ-I.I)=BC(J.I) 
1570 


683 
BBCJ-l.2)=BC(J-I.l) 
1571 
684 
IF C J - 2) 
900. 900.916 
1572 
9 0 0 
I F( BC C 1 • 1 ) -BC C 2. 1 )) 
91 5. 916 • 91 6 
1 57 3 
915 
B8(1.2)=8C(1.I) 
1574 
916 
T(J)=SAYEI4 
1575 


CWUJ.l )=SAYEI5 
1576 


C 
CALCULATE TEMPORARY FLUXS 
1577 
831 
JSTR=J+2 
1578 
832 
IFeJSTR-LZR) 
834.834.833 
1579 
833 
JSTR=LZR 
1580 


834 
DO 
842 N=I.NSZ 
1581 


159 


835 
FITCJSTR+I.N)=FICJSTR+l.N) 
1582 
836 
00 838 
L=I.JSTR 
1583 
837 
K=JSTR-L+I 
1584 
838 
F1TCK.N)=FITCK+I.N). ZICK+I.N)+BBCK.I).C ACK+I.N»-2.0. SCK.N) 
1585 
838 
2. ~CK+l.N) 
:586 


839 
FOTCI.N)=FITCI.N). ZICI.N)+8BCI.2).C ACI.N))+2.0. Sel.N) 
1587 
839 2. ~CI.N) 
1588 


840 
00 841 
K=2.JSTR 
1589 
841 
FOTCK.N)=FOTCK-I.N).CR2CK-I)/R2CK)). lOCK.N) 
1590 


841 
2+BBCK.2).C ACK.N»+2.0. SCK.N). ~CK.N) 
1591 
841 
3+CI.0-CR2CK-I)/R2eK»).FITCK.N)* lICK.N) 
:592 


842 
CONTINUE 
:593 


C 
SPECTRAL SUMMATION 
1594 
726 
FOTJP2=0.0 
1595 
727 
FOTJPl=O.O 
1596 
728 
FOTJl=O.O 
1597 
729 
FOTJMI=O.O 
1598 


730 
FITJZ=O.O 
1599 
731 
FITJMI=O.O 
1600 


732 
FITJM2=0.O 
1601 


733 
FITJM3=0.O 
1602 


734 
00 742 
N=I.NSZ 
1603 


735 
FOTJP2=FOTJP2 + FOTCJ+2.N) , 
1604 


736 
FOTJPI=FOTJPI + FOTCJ+I.N) 
1605 


737 
FOTJl=FOTJZ + FOTCJ.N) 
1606 


738 
FOTJMI=FOTJMI + FOTeJ-I.N) 
1607 


739 
FITJZ=FITJZ + FITCJ.N) 
1608 
740 
FITJMl=FITJMI + FITCJ-I.N) 
1609 


741 
FITJM2=FITJM2 + FITCJ-2.N) 
1610 


742 
FITJM3=FITJM3 + FITeJ-3.N) 
1611 


FOCJ.2)=FOTJl 
1612 


FleJ.2)=FITJZ 
1613 


C 
CALCULATE FLUX DERIVATIVES 
1614 
743 
IFeJ-LZR+l) 
744.745.746 
1615 
744 
FMMCJ+2)=eFOTJP2 -FOSeJ+2»)/eTPCJ)-T(J» 
1616 
745 
FMeJ+l)=CFOTJPl -FOSeJ+l»/(TPCJ)-TeJ» 
1617 
746 
FZeJ)=eFOTJZ-FOSCJ»/eTPeJ)-TCJ» 
1618 
747 
HZeJ)=eFITJZ-FISeJ»/eTPeJ)-TeJ» 
1619 
748 
IFeJ-l) 
755.755.749 
1620 


749 
FPeJ-I)=(FOTJMl -FOSeJ-I»/(TPeJ)-TeJ» 
1621 


750 
HP(J-I)=eFITJMI-FISCJ-l»/CTP(J)-TeJ» 
1622 
751 
IFeJ-2) 
755.755.752 
1623 
752 
HPpeJ-2)=(FITJM2-FIS(J-2»/eTPeJ)-T(J» 
1624 
753 
IFeJ-3) 
755.755.754 
1625 
754 
HPPPCJ-3)=eFITJM3-FISeJ-3»/CTPeJ)-TeJ» 
1626 


C 
RETURN RIPPLE ZONE PARAMETERS TO NORMAL VALUES 
1627 
755 
00 768 
N=I.NSZ 
1628 


BBeJ.N)=SAVEI 
1629 


BBeJ-I.N)=SAVE2 
1630 


SeJ.N)=SAVE3 
1631 
SeJ-l.N)=SAVE4 
1632 
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OTAUCJ.Nl=SAYE5 
1633 


ACJ.Nl=SAYE6 
1634 


VCJ.Nl=SAYE7 
1635 


ZICJ.Nl=SAYE9 
1636 


ZJ(J-l.Nl=SAYE9 
1637 


ZOeJ.Nl=SAYE10 
1639 


lOeJ+I.Nl=SAYEII 
1639 


9CeJ.Nl=SAYEI2 
1640 
TSeJl=SAYEI3 
1641 


ZOCJ+2.Nl=SAYEI6 
1642 
BBeJ.2)=SAYE17 
1643 


BBeJ-I.2)=SA.YEI9 
1644 
769 
CONTINUE 
1645 
769 
CONTINUE 
1646 


C GENERALIZED COEFFICIENT ROUTINE 
DEC 6.1963 
1647 


770 
DO 799 
K=I,LZR 
1648 
771 
IFCK-I) 
772.772.775 
1649 
772 
OMMe 1 )=0. 0 
1650 
01=0.0 
1651 
773 
OM(I)=O.O 
1652 
774 
GO TO 780 
1653 
775 
01=1.0 
1654 


IFCK-2) 
777.776.779 
1655 
776 
OMM(2)=0.0 
1656 
777 
GO TO 779 
1657 
778 
OMMCK)=CO.5*R2CK)/ZMASCK»*FMMCK) 
1658 
779 
1-01*CO.5*R2CK-I)/ZMASCK»*~MCK-I) 
1659 
779 
OMCK)= CO.5*R2CK)/ZHASCK»wFHCK) 
1660 


779 
1-01*CO.5*R2CK-11/ZHASCK»*CFZCK-I)-HZCK-I» 
1661 


790 
OZCK)= CO.5*R2CK)/ZMASCK»*CFZCK1-HZCK» 
1662 


790 
1-01*CO.5*R2CK-l1/ZMASCK»*£FPCK-ll-HPCK-I» 
1663 


780 
2+COETCK)/OT)+COPTCK)/C2.*OT»*CVCK)-VRCK» 
1664 
781 
OPCK)= CO.5*R2CK)/ZMASCK»*CFPCK)-HP(K» 
1665 
781 
1-01*CO.5*R2(K-l)/ZHASCK»*C-HPPCK» 
1666 
782 
OPPCK)=CO.5*R2CK)/ZHASCK»*C-HPPCK» 
1667 
782 
1-01*CO.5*R2CK-l)/ZMASCK»*C-HPPPCK-l» 
1668 
799 
CONT INUE 
1669 


790 
OPPCLZR-l )=0. 0 
1670 


791 
OPCLZR)=O.O 
1671 
792 
OPPCLZR)=O.O 
1672 


C BLOCK 900 SOLUTION OF HATRIX 
1673 


900 
CN(1)=DPC1)/DZCl) 
1674 


901 
ON(I)=DPP(I)/OZ(I) 
1675 


902 
EN(ll=-RESOUECI)/DZCIl 
1676 


903 
BN(2)=OZC21-DHC21*CNCI) 
1677 


904 
CNC21=COP(2)-DM(2)tONCI)/BNC2) 
1679 


905 
DN(2)=OPPC2)/BN(2) 
1679 


906 
EN(2)=-CRESDUEC21+0MC21tENCI1)/BNC2) 
1690 
807 
DO 811 
K=3.LZR 
1681 


808 
BNCK)=eozeKl-OMCKltCNCK-Il-OMMCK)tDNeK-21+0MMCK1*CNCK-21*CNCK-I» 
1682 


809 
CNCK)=COPCKl-OMCKltONCK-I)+OMMCK)*CNCK-2)*ONCK-I»/BNCK) 
1683 
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810 
DNCK)=DPPCK)/8NCK) 


811 
ENCK)=-CRESDUECK)+DMCK).ENCK-l)-DMMCK).CNCK-2).ENCK-I)+DMMCK).ENCK 


811 
1-2»/BNCK) 


812 
QNCLZR-I)=O.O 


813 
BNCLZR)=DMCLZR)-DMMCLZR).CNCLZR-2) 


814 
CNCLZR)=CDZCLZR)-DMMCLZR).DNCLZR-2»/BNCLZR) 


815 
DNCLZR)=O.O~ 


816 
ENCLZR)=-CRESDUECLZR)+DMMCLZR).ENCLZR-2)/BNeLZR) 


819 
RETURN 


820 
END 


C 
VERSION C ENTRY ROUTINE 
STARTING MODEL IS READ IN FROM DATA CARDS GENERATED 
BY THE X-RAY DEPOSIT ROUTINE ON UNIYAC 1107 
• 
FORTRAN 
SUBROUTINE ENTRY 
C 
STANDARD DIMENSION AND COMMON STATEMENTS 
AS IN MAIN PROGRAM 


2 
DCI)=O.O 
TEE=3.0E+04 
READ INPUT TAPE 5. 9277.RUNIDI.RUNID2.TReNNI).TRCNN2).ENERGY.TIME. 
IFX.CK.RCK).DRCK).TCK).RHOZCK).ueK).ECK).RUNIOl.RUNI02.K=I.IOO) 


9277 
FORMATC2A6.IP5EI3.3/CI4.IP6EII.3.2A5» 
LZ=IOO 
YD=CS 
WKT=ENERGY/4.18E+19 
TOO2=7100.0.EXPFceO.43429.LOGFCRHOZ(99)/1.293E-03»/5.3) 
TON2=15000.0.EXPF(CO.43429.LOGFeRHOZe99)/1.293E-03»/4.8) 
,~ 
1019 00 1033 K=I.I00 
1026 VCK)=I.O/RHOZCK) 
YReK)=VCK) 
1027 L=K 
TSCK)=T(K) 
R2CL)=ReUu2 
1029 RzeL)=Reu 
1030 RReu=1. 0 
TRCK)=T(K) 
ZMASCK)=CRHOZeK)/3.0).ceRZeL) •• 3)-CRZeL-l) •• 3» 
OU=UCU-UCL-I) 
QCK)=-CVD/C2.0.YCK»).OU.ABSFCOU) 
Q(I)=-CVD/C2.0.YCI»).UCI).ABSFCU(I» 
IFCQCK» 
2233.1033.1033 


2233 
QCtO =0.0 
1033 
CONT INUE 
ZMASCI)=RHOZ(I).CRCI)·.3)/3.0 
1034 WRITE OUTPUT TAPE 6.1035.TDN2.TD02.RN2.R02.WKT 
.CK.RCK).RHOZCK). 
1034 lVCKJ.TCK).UCK).QCK).ZMASCK).DRCK).DPPCK).DTAUCK.l).TAUCK.l).K=I. 
1100) 
1035 FORMATCIH .10HINPUT DATA. IOX.5HTDN2=OPF7. 1. 10X.5HTD02=OPF6. 1.5X. 


1035 14HRN2=IPE9.2.5X.4HR02=IPE9.2.5X.4HWKT=OPF7.2114H 
K .2X.4HRCL) 
I 
.7X.7HRHOZCK).4X.4HYCK).7X.4HTCK1.7X.4HUCL).7X.4HQCK1.7X.7HZMAS 
ICK).4X.5HDRCK).7X.6HENGCK).5X.7HDTAUCK).4X.3HTAUIICI4.IPIIEII.4» 
T( 101 )=TC I 00) 
1000 
RETURN 
END 
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